Tutorial: Data-driven emission quantification of hot spots¶

This notebook demonstrates how to quantify CO\(_2\) and NO\(_x\) emissions from point sources using synthetic CO2M observations for a power plant and for a city. The data files used in this tutorial are a subset of the SMARTCARB dataset and can be found in the ddeq.DATA_PATH folder. The full SMARTCARB dataset is avalable here: https://doi.org/10.5281/zenodo.4048227

[1]:

import os
import ucat

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import pandas as pd
import xarray as xr

from ddeq import DATA_PATH
from ddeq.smartcarb import DOMAIN
import ddeq

# Get optimal coordinate reference system for computing Easting and Northing in DOMAIN
CRS = ddeq.misc.get_opt_crs(DOMAIN)

2025-12-04 09:14:44.871780065 [W:onnxruntime:Default, device_discovery.cc:164 DiscoverDevicesForPlatform] GPU device discovery failed: device_discovery.cc:89 ReadFileContents Failed to open file: "/sys/class/drm/card0/device/vendor"

Read list of point sources in the SMARTCARB model domain from “sources-smartcarb.csv” file. The format of the xr.Dataset is used by plume detection and emission quantification code internally to identify the point sources.

The dataset includes names and locations of sources as well as annual mean CO\(_2\) and NO\(_x\) emissions used in the SMARTCARB simulations. Note that the true emissions in the simulations varying temporally. The standard deviation of the temporal variability is given as precision.

[2]:

# list of point sources
sources = ddeq.sources.read_smartcarb()
sources

Plume detection algorithm¶

Plumes are regions how satellite pixels where CO2/NO2 values are significantly enhanced above the background \begin{equation} SNR = \frac{X - X_\mathrm{bg}}{\sqrt{\sigma_\mathrm{random}^2 + \sigma_\mathrm{sys}^2}} \geq z_\mathrm{thr} \end{equation} The value \(X\) is computed by applying a Gaussian filter (other filters are possible) with size filter_size (default: 0.5 pixels). The background \(X_\mathrm{bg}\) is computed using a median filter (size = 100 pixels). The threshold \(z_\mathrm{thr}\) is computed for z-statistics using a probability \(q\) (default: 0.99). Pixels for which above equation is true, are connected to regions using a labeling algorithm considering (horizontal, vertical and diagonal neighbors). Regions that overlap that are within the radius defined in sources of a point sources are assigned to the source. A region can be assigned to more than one source (overlapping plumes).

[8]:

filename = os.path.join(DATA_PATH, 'Sentinel_7_CO2_2015042311_o1670_l0483.nc')

data = ddeq.smartcarb.read_level2(filename, co2_noise_scenario='low',
                                  no2_noise_scenario='high',
                                  co_noise_scenario='low',
                                  only_observations=False)

[9]:

data = ddeq.dplume.detect_plumes(data, sources, variable='NO2', variable_std='NO2_std',
                                 filter_type='gaussian', filter_size=0.5, crs=CRS)

The code computes several new fields that are added to the provided data dataset. The detected plumes are stored in the detected_plume data array (dims: nobs, nrows, source) where the length of source is equal to the number of detected plumes.

The plume detection can be visualized using the ddeq.vis.show_level2 function:

[10]:

fig = ddeq.vis.show_level2(data, data['NO2'], gas='NO2', domain=DOMAIN, sources=sources,
                           winds=winds, do_zoom=False, show_clouds=True)

_images/tutorial-introduction-to-ddeq_24_0.png

Exercise¶

Detect emission plumes Berlin and Jänschwalde using low-noise CO2 observations.
Increase the size of the Gaussian filter to increase number of detected pixels.
Visualize the results using ddeq.vis.show_level2

[ ]:

Center lines and polygons¶

To estimate emissions for detected plumes, the following code fits a center curve for each detected plume. The code also adds across- and along-plume coordinates (xp and yp). The computation of plume coordinates can result in multiple solutions when the center line is strongly curved or the plume is small.

[11]:

# Read SMARTCARB Level-2 file
filename = os.path.join(DATA_PATH, 'Sentinel_7_CO2_2015042311_o1670_l0483.nc')
data = ddeq.smartcarb.read_level2(
    filename, co2_noise_scenario='low', no2_noise_scenario='high',
    co_noise_scenario='low', only_observations=False
)

# Detect plumes using NO2 observations
data = ddeq.dplume.detect_plumes(
    data, sources.sel(source=['Berlin', 'Janschwalde']), crs=CRS,
    variable='NO2', variable_std='NO2_std',
    filter_type='gaussian', filter_size=0.5
)

# Fit center curve to detected plumes
data = ddeq.curves.fit_to_detections(data, n_nodes=3, force_origin=True, use_weights=True)

# Compute natural coordinates of along and across plume direction for each plume:
data = ddeq.curves.compute_natural_coords(data)

# Define areas around plume
data = ddeq.curves.compute_plume_areas(data)

The following code shows the result:

[12]:

ddeq.vis.show_level2(
    data, 'NO2', gas="NO2", domain=DOMAIN, winds=winds, do_zoom=True,
    show_clouds=True, draw_gridlines=True, crs=CRS
);

_images/tutorial-introduction-to-ddeq_31_0.png

Prepare emission quantification¶

To prepare for estimating the emissions the following code computes the CO2 and NO2 background field, the plume signals and converts to mass columns in kg/m² using the ucat Python package.

[13]:

for gas in ["CO2", "NO2"]:

    # estimate background
    data = ddeq.background.estimate(data, gas)

    # compute CO2/NO2 enhancement
    data = ddeq.emissions.compute_plume_signal(data, gas)

    # convert ppm to kg/m2
    for variable in [
        gas,
        f"{gas}_estimated_background",
        f"{gas}_minus_estimated_background",
    ]:
        ddeq.emissions.convert_units(data, gas, variable)

It is possible to visualize different variables using the variable parameter:

[14]:

ddeq.vis.show_level2(
    data, 'CO2_minus_estimated_background_mass', gas='CO2', domain=DOMAIN,
    winds=None, do_zoom=True, show_clouds=True, draw_gridlines=True,
    vmin=-20e-3, vmax=40e-3,
    label='CO$_2$ enhancement [kg m$^{-2}$]', crs=CRS
);

_images/tutorial-introduction-to-ddeq_36_0.png

Integrated mass enhancement¶

The following code uses integrated mass enhancement for computing CO\(_2\) emissions. First, the wind field is extracted from the SMARTCARB dataset at each source location. Second, the IME method is applied to estimate the emissions.

[20]:

winds = ddeq.wind.read_smartcarb(time, sources.lon, sources.lat, data_path=DATA_PATH)

[21]:

results = ddeq.ime.estimate_emissions(data, winds, sources, gas='CO2')

print(' ' * 10, '\tEstimate\tTrue')

for name, source in sources.groupby('source'):
    if name in data.source:
        Q = results['CO2_emissions'].sel(source=name).values
        Q_true = ddeq.smartcarb.read_true_emissions(
            time=pd.Timestamp(data.time.values), gas='CO2', source=name
        ).mean()
        print(
            f'{name:10s}\t'
            f'{ucat.convert_mass_per_time_unit(Q, "kg/s", "Mt/a"):.1f} Mt/a\t'
            f'{ucat.convert_mass_per_time_unit(Q_true, "kg/s", "Mt/a"):.1f} Mt/a'
        )

                Estimate        True
Berlin          19.3 Mt/a       23.4 Mt/a
Janschwalde     43.4 Mt/a       42.4 Mt/a

The results can be visualized with the following function:

[22]:

Q_true = ddeq.smartcarb.read_true_emissions(
    time=pd.Timestamp(data.time.values),
    gas='CO2',
    source='Janschwalde'
).mean()

with xr.set_options(keep_attrs=True):
    fig = ddeq.vis.plot_ime_result(
        'CO2',
        data, winds, results,
        source='Janschwalde',
        domain=DOMAIN, crs=CRS,
        true_emissions=Q_true,
        do_zoom=True,
    )

_images/tutorial-introduction-to-ddeq_55_0.png

[23]:

results = ddeq.ime.estimate_emissions(data, winds, sources,
                                      gas='NO2', decay_time=4*60**2)

results = ddeq.emissions.convert_NO2_to_NOx_emissions(results, f=1.32)

print(' ' * 10, '\tEstimate\tTrue')

for name, source in sources.groupby('source'):
    if name in data.source:
        Q = results['NOx_emissions'].sel(source=name).values
        Q_true = ddeq.smartcarb.read_true_emissions(
            time=pd.Timestamp(data.time.values), gas='NO2', source=name
        ).mean()
        print(
            f'{name:10s}\t'
            f'{ucat.convert_mass_per_time_unit(Q, "kg/s", "kt/a"):.1f} Mt/a\t'
            f'{ucat.convert_mass_per_time_unit(Q_true, "kg/s", "kt/a"):.1f} Mt/a'
        )

                Estimate        True
Berlin          31.2 Mt/a       25.2 Mt/a
Janschwalde     46.1 Mt/a       34.2 Mt/a

[24]:

Q_true = ddeq.smartcarb.read_true_emissions(
    time=pd.Timestamp(data.time.values),
    gas='NOx',
    source='Janschwalde'
).mean()

with xr.set_options(keep_attrs=True):
    fig = ddeq.vis.plot_ime_result(
        'NO2',
        data, winds, results,
        source='Janschwalde',
        domain=DOMAIN, crs=CRS,
        true_emissions=Q_true,
        do_zoom=True,
    )

_images/tutorial-introduction-to-ddeq_57_0.png

[ ]:

Tutorial: Data-driven emission quantification of hot spots¶

Synthetic satellite observations¶

Exercise¶

Wind fields¶

Plume detection algorithm¶

Exercise¶

Center lines and polygons¶

Prepare emission quantification¶

Cross-sectional flux method¶

True emissions¶

Cross-sectional flux method for point source¶

Cross-sectional flux method for a city¶

Exercise¶

Integrated mass enhancement¶