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.
[1]:
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()
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 = ccrs.epsg(22293)
[2]:
data = xr.open_dataset(filename)
[3]:
fig = ddeq.vis.visualize(
data, 1e6 * data['NO2'],
sources=sources,
vmin=0,
vmax=300,
label="NO$_2$ vertical column density [µmol m$^{-2}$]",
units='umol m-2'
)

Cross-sectional flux method#
[4]:
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)
data, curves = ddeq.plume_coords.compute_plume_line_and_coords(
data, crs=CRS, radius=25e3, plume_area='area'
)
data = ddeq.emissions.prepare_data(data, 'NO2')
[5]:
winds = ddeq.wind.read_at_sources(data.time,
sources.sel(source=['Matimba']),
timesteps=12,
era5_prefix='Matimba',
data_path=ddeq.DATA_PATH)
[6]:
results_csf = ddeq.csf.estimate_emissions(data, winds, sources, curves, 'NO2',
f_model=2.24, crs=CRS)
[7]:
fig = ddeq.vis.show_detected_plumes(
data, curves, 1e6 * data['NO2'], gas='NO2', ld=results_csf,
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,
);

[8]:
fig, ax = plt.subplots(1, figsize=(8,2))
ddeq.vis.plot_along_plume(ax, 'NO2', results_csf['Matimba'])
ax.set_xlim(right=207)
ax.set_ylim(0,150)
plt.tight_layout()
ax.legend().remove()
ax.legend(loc=1, ncol=3)
[8]:
<matplotlib.legend.Legend at 0x7f11de482050>

[9]:
fig, axes = plt.subplots(1,3, figsize=(8,2), sharey=True)
for i, ax in zip([1,9,19], axes):
r = results_csf['Matimba'].isel(along=i)
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()

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

Gaussian plume inversion#
[11]:
priors = {'Matimba': {
'NO2': {'Q': 3.0, # kg/s
'tau': 4*60**2 # seconds
}
}}
data, results_gauss = ddeq.gauss.estimate_emissions(data, winds, sources, curves,
['NO2'], priors=priors,
fit_decay_times=True)
[12]:
results_gauss = ddeq.emissions.convert_NO2_to_NOx_emissions(results_gauss, f=2.38)
[13]:
fig = ddeq.vis.plot_gauss_result(data, results_gauss, ['Matimba'],
'NO2', curves, crs=CRS, vmin=0, vmax=10e-6)

[14]:
data['NO2_plume_model'] = xr.DataArray(
ucat.convert_columns(data['NO2_plume_model_mass'], 'kg m-2', 'mol m-2',
molar_mass='NO2'),
dims=data.NO2_plume_model_mass.dims
)
[15]:
fig = ddeq.vis.show_detected_plumes(
data, curves,
1e6 * data['NO2_plume_model'].sel(source='Matimba'),
gas='NO2',
vmin=0,
vmax=300,
label="NO$_2$ column [µmol m$^{-2}$]",
winds=winds,
domain=DOMAIN,
do_zoom=False, show_clouds=False, crs=CRS,
figwidth=4, add_polygons=False
);
xp, yp = results_gauss[f'NO2_curve'].sel(source='Matimba').T
x, y = ddeq.misc.transform_coords(xp, yp, CRS, ccrs.PlateCarree(), use_xarray=False)
fig.axes[0].plot(x, y, 'r', transform=ccrs.PlateCarree())
[15]:
[<matplotlib.lines.Line2D at 0x7f11e1a8a740>]
