Example: TROPOMI NO\(_2\) plume

An example demonstrating the application of the CSF and Gaussian plume inversion for estimating NOx emissions from the Matimba/Medupi power plant using TROPOMI NO2 images.

import os

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

import ddeq

sources = ddeq.misc.read_point_sources()
sources = sources.sel(source=['Matimba'])

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

DOMAIN = ddeq.misc.Domain('Matimba', 25, -25.2, 29, -22.9)
CRS = ddeq.misc.get_opt_crs(DOMAIN)

CURVE_STYLE = "BEZIER"     # "POLY" or "BEZIER"
DETECTION = "THR"         # "THR" or "WIND"
BACKGROUND = "LINEAR"      # "LINEAR" or "SMOOTH"

2025-12-04 09:16:01.001243703 [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"

Sentinel-5P/TROPMI NO2 dataset

data = xr.open_dataset(filename)

fig = ddeq.vis.show_level2(
    data, 1e6 * data['NO2'],
    sources=sources,
    vmin=0,
    vmax=300,
    label="NO$_2$ vertical column density [µmol m$^{-2}$]", do_zoom=True,
    units='umol m-2'
)

ERA5 wind files

The code below would download ERA5 single level and pressure level. Since the files are already in ddeq.DATA_PATH, no files will be downloaded. If you want to download own ERA5 data, please the change data_path to the location on your system, where you like to save your files.

sng_filename = ddeq.era5dl.download_single_lvl(
    time=pd.Timestamp(data.time.values),
    area=DOMAIN.extent,
    timesteps=24,
    data_path=ddeq.DATA_PATH,
    prefix="Matimba"
)

lvl_filename = ddeq.era5dl.download_pressure_lvl(
    time=pd.Timestamp(data.time.values),
    area=DOMAIN.extent,
    timesteps=24,
    data_path=ddeq.DATA_PATH,
    prefix="Matimba"
)

sng_filename, lvl_filename

('/home/docs/checkouts/readthedocs.org/user_builds/ddeq/envs/stable/lib/python3.12/site-packages/ddeq/data/Matimba_ERA5-sl-20210725.nc',
 '/home/docs/checkouts/readthedocs.org/user_builds/ddeq/envs/stable/lib/python3.12/site-packages/ddeq/data/Matimba_ERA5-pl-20210725.nc')

Cross-sectional flux method

# Effective wind at source locations using GNFR-A weighted winds:
winds = ddeq.era5.read(
    sng_filename,
    lvl_filename,
    method="gnfra",
    sources=sources,
    times=data.time
)
winds

<xarray.Dataset> Size: 96B
Dimensions:          (time: 1, source: 1)
Coordinates:
  * time             (time) datetime64[ns] 8B 2021-07-25T11:44:52.595066640
  * source           (source) object 8B 'Matimba'
    number           int64 8B 0
    expver           (time) <U4 16B '0001'
    lon              (source) float64 8B 27.61
    lat              (source) float64 8B -23.67
Data variables:
    U                (time, source) float64 8B -5.856
    V                (time, source) float64 8B -2.344
    speed            (time, source) float64 8B 6.308
    speed_precision  (time, source) float64 8B 1.0
    direction        (time, source) float64 8B 68.18
Attributes:
    CREATOR:       ddeq.era5
    DATE_CREATED:  2025-12-04 09:16
    ORIGIN:        ERA-5
    METHOD:        gnfra
    DESCRIPTION:   Weighted average for U- and V-wind using GNFR-A emission p...

xarray.Dataset

• Dimensions:
  • time: 1
  • source: 1
• Coordinates: (6)
  • time
    (time)
    datetime64[ns]
    2021-07-25T11:44:52.595066640

    array(['2021-07-25T11:44:52.595066640'], dtype='datetime64[ns]')

  • source
    (source)
    object
    'Matimba'

    array(['Matimba'], dtype=object)

  • number
    ()
    int64
    0

    long_name :

    ensemble member numerical id

    units :

    1

    standard_name :

    realization

    array(0)

  • expver
    (time)
    <U4
    '0001'

    array(['0001'], dtype='<U4')

  • lon
    (source)
    float64
    27.61

    name :

    longitude of point source

    array([27.610556])

  • lat
    (source)
    float64
    -23.67

    name :

    latitude of point source

    array([-23.668333])

• Data variables: (5)
  • U
    (time, source)
    float64
    -5.856

    GRIB_paramId :

    131

    GRIB_dataType :

    an

    GRIB_numberOfPoints :

    170

    GRIB_typeOfLevel :

    isobaricInhPa

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_uvRelativeToGrid :

    0

    GRIB_NV :

    0

    GRIB_Nx :

    17

    GRIB_Ny :

    10

    GRIB_cfName :

    eastward_wind

    GRIB_cfVarName :

    u

    GRIB_gridDefinitionDescription :

    Latitude/Longitude Grid

    GRIB_iDirectionIncrementInDegrees :

    0.25

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.25

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    -22.95

    GRIB_latitudeOfLastGridPointInDegrees :

    -25.2

    GRIB_longitudeOfFirstGridPointInDegrees :

    25.0

    GRIB_longitudeOfLastGridPointInDegrees :

    29.0

    GRIB_missingValue :

    3.4028234663852886e+38

    GRIB_name :

    U component of wind

    GRIB_shortName :

    u

    GRIB_totalNumber :

    0

    GRIB_units :

    m s**-1

    long_name :

    U component of wind

    units :

    m s-1

    standard_name :

    eastward_wind

    description :

    GNFR-A weighted u-wind

    array([[-5.85571034]])

  • V
    (time, source)
    float64
    -2.344

    GRIB_paramId :

    132

    GRIB_dataType :

    an

    GRIB_numberOfPoints :

    170

    GRIB_typeOfLevel :

    isobaricInhPa

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_uvRelativeToGrid :

    0

    GRIB_NV :

    0

    GRIB_Nx :

    17

    GRIB_Ny :

    10

    GRIB_cfName :

    northward_wind

    GRIB_cfVarName :

    v

    GRIB_gridDefinitionDescription :

    Latitude/Longitude Grid

    GRIB_iDirectionIncrementInDegrees :

    0.25

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.25

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    -22.95

    GRIB_latitudeOfLastGridPointInDegrees :

    -25.2

    GRIB_longitudeOfFirstGridPointInDegrees :

    25.0

    GRIB_longitudeOfLastGridPointInDegrees :

    29.0

    GRIB_missingValue :

    3.4028234663852886e+38

    GRIB_name :

    V component of wind

    GRIB_shortName :

    v

    GRIB_totalNumber :

    0

    GRIB_units :

    m s**-1

    long_name :

    V component of wind

    units :

    m s-1

    standard_name :

    northward_wind

    description :

    GNFR-A weighted v-wind

    array([[-2.34447273]])

  • speed
    (time, source)
    float64
    6.308

    units :

    m s-1

    array([[6.3076062]])

  • speed_precision
    (time, source)
    float64
    1.0

    units :

    m s-1

    array([[1.]])

  • direction
    (time, source)
    float64
    68.18

    units :

    °

    array([[68.18013209]])

• Attributes: (5)

  CREATOR :

  ddeq.era5

  DATE_CREATED :

  2025-12-04 09:16

  ORIGIN :

  ERA-5

  METHOD :

  gnfra

  DESCRIPTION :

  Weighted average for U- and V-wind using GNFR-A emission profile.

if DETECTION == "THR":

    var_sys = (ucat.convert_columns(0.5e15, 'cm-2', 'mol m-2',
                                    molar_mass='NO2'))**2

    data = ddeq.dplume.detect_plumes(data, sources,
                                     variable='NO2',
                                     variable_std='NO2_std',
                                     var_sys=var_sys,
                                     filter_type='gaussian',
                                     filter_size=0.5, crs=CRS)
elif DETECTION == "WIND":
    data = ddeq.dplume.detect_from_wind(data, sources, winds, dmax=30e3, crs=CRS)
else:
    raise ValueError

fig = ddeq.vis.show_level2(
    data, 1e6 * data['NO2'], gas='NO2',
    vmin=0,
    vmax=300,
    label="NO$_2$ columns [µmol m$^{-2}$]",
    winds=winds,
    domain=DOMAIN,
    do_zoom=False, show_clouds=False, crs=CRS,
    figwidth=4,
);

variable = "NO2_mass"
background = "linear"
ddeq.emissions.convert_units(data, "NO2", "NO2")

if DETECTION == "THR":
    data = ddeq.curves.fit_to_detections(data, n_nodes=3, force_origin=True, use_weights=True)

data = ddeq.curves.compute_natural_coords(data)
data = ddeq.curves.compute_plume_areas(data, plume_width=120e3)

# estimate NO2 enhancement
ddeq.background.estimate(data, "NO2")
ddeq.emissions.compute_plume_signal(data, "NO2")

ddeq.emissions.convert_units(data, "NO2", "NO2_estimated_background")
ddeq.emissions.convert_units(data, "NO2", "NO2_minus_estimated_background")

dx = None if DETECTION == "WIND" else 10e3

results_csf = ddeq.csf.estimate_emissions(data, winds, sources, 'NO2',
                                          xmin=10e3, dx=dx,
                                          f_model=2.24, crs=CRS,
                                          variable=variable, background=background
                                         )

/home/docs/checkouts/readthedocs.org/user_builds/ddeq/envs/stable/lib/python3.12/site-packages/ddeq/csf.py:1098: FutureWarning: In a future version of xarray the default value for data_vars will change from data_vars='all' to data_vars=None. This is likely to lead to different results when multiple datasets have matching variables with overlapping values. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set data_vars explicitly.
  polygons = xr.concat(

ddeq.vis.plot_csf_result(['NO2'], data, winds, results_csf, 'Matimba',
                         vmins=[0], vmaxs=[200e-6], crs=CRS);

fig = ddeq.vis.show_level2(
    data,
    1e6 * data['NO2'],
    gas='NO2',
    vmin=0,
    vmax=300,
    label="NO$_2$ columns [µmol m$^{-2}$]",
    winds=winds,
    domain=DOMAIN,
    do_zoom=False,
    show_clouds=False,
    crs=CRS,
    figwidth=4,
);

fig, ax  = plt.subplots(1, figsize=(8,2))
ddeq.vis.plot_along_plume(ax, 'NO2', results_csf.sel(source='Matimba'))
ax.set_xlim(right=207)
ax.set_ylim(0,150)
plt.tight_layout()

ax.legend().remove()
ax.legend(loc=1, ncol=3)

<matplotlib.legend.Legend at 0x7d8b99ba1520>

fig, axes  = plt.subplots(1,3, figsize=(8,2), sharey=True)


for i, ax in zip([1,9,-1], axes):
    try:
        r = results_csf.sel(source='Matimba').isel(polygon=i)
    except IndexError:
        continue
    ddeq.vis.plot_across_section(r, gases=['NO2'], method='gauss', ax=ax,
                                 legend='simple')

    ax.text(-36, 19.5, '%d-%d km' % (r.xa/1e3, r.xb/1e3), ha='left', va='top')

for ax in axes[1:]:
    ax.set_ylim(-1, 20)
    ax.set_ylabel('')


plt.tight_layout()

ddeq.vis.plot_csf_result(['NO2'], data, winds, results_csf, 'Matimba',
                         vmins=[0], vmaxs=[200e-6], crs=CRS);

Gaussian plume inversion

priors = {'Matimba': {
    'NO2': {'Q': 3.0,       # kg/s
            'tau': 4*60**2  # seconds
           }
}}
data, results_gauss = ddeq.gauss.estimate_emissions(data, winds, sources,
                                                    ['NO2'], priors=priors,
                                                    fit_decay_times=True)

results_gauss = ddeq.emissions.convert_NO2_to_NOx_emissions(results_gauss, f=2.38)

fig = ddeq.vis.plot_gauss_result(data, results_gauss, ['Matimba'],
                                 'NO2', crs=CRS, vmin=0, vmax=10e-6)

IME method

results_ime = ddeq.ime.estimate_emissions(
    data,
    winds,
    sources,
    "NO2",
    decay_time=4.0*60**2
)
results_ime = ddeq.emissions.convert_NO2_to_NOx_emissions(results_ime, f=1.32)

ddeq.vis.plot_ime_result(
    "NO2",
    data,
    winds,
    results_ime,
    "Matimba",
    crs=CRS,
    vmin=0,
    vmax=300e-6,
);

LCSF method

Quick example applying the LCSF method to TROPOMI NO2 data

# Convert from mol/m2 to molec/cm2 for LCSF method
data["NO2"] = ucat.convert_columns(data["NO2"], "mol/m2", "cm-2")
data["NO2_std"] = ucat.convert_columns(data["NO2_std"], "mol/m2", "cm-2")

wind_field = ddeq.era5.read(
    sng_filename,
    lvl_filename,
    method="gnfra",
    times=data.time
)

lcs_params = {}
lcs_params['use_prior'] = True
lcs_params["verbose"] = True

lcs_params['n_min_fit_pts'] = 10
lcs_params['fit_pt_slice_width'] = 15

lcs_params['f_NOx_NO2'] = 3.5
lcs_params['NOx_NO2_tau_depletion'] = 4.0 * 60**2

lcsf_results = ddeq.lcsf.estimate_emissions(data, wind_field, sources, ['NO2'], priors=priors, lcs_params=lcs_params, all_diags=True)
print(f"NOx emissions: {np.nanmedian(lcsf_results['NO2_emissions']):.1f} kt NO2 / a")

LCSF parameters
{'use_prior': True, 'verbose': True, 'n_min_fit_pts': 10, 'fit_pt_slice_width': 15, 'f_NOx_NO2': 3.5, 'NOx_NO2_tau_depletion': 14400.0, 'min_estim_emis': None, 'max_estim_emis': None, 'Enhancement_thresh': 1.0, 'max_sigma_Gauss': 5.0, 'max_rel_std_emis': 1.0, 'window_length': 100}

**************

Estimating line densities and emissions for source Matimba

10 LD found !
NO2 median emission estimate  101.68162
NOx emissions: 101.7 kt NO2 / a

fig = ddeq.vis.plot_lcsf_result("Matimba", lcsf_results, data, sources,
                                gases=['NO2'])

/home/docs/checkouts/readthedocs.org/user_builds/ddeq/envs/stable/lib/python3.12/site-packages/ddeq/vis.py:1097: UserWarning: To plot the polygons, a coordinate reference system (crs) must be provided.
  warnings.warn(