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')
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[4], line 15
4 data = ddeq.dplume.detect_plumes(data, sources,
5 variable='NO2',
6 variable_std='NO2_std',
7 var_sys=var_sys,
8 filter_type='gaussian',
9 filter_size=0.5)
11 data, curves = ddeq.plume_coords.compute_plume_line_and_coords(
12 data, crs=CRS, radius=25e3, plume_area='area'
13 )
---> 15 data = ddeq.emissions.prepare_data(data, 'NO2')
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/emissions.py:74, in prepare_data(data, gas)
70 molar_mass = str(gas)
72 psurf = data['psurf'] if 'psurf' in data else data['suface_pressure']
73 data[name] = xarray.DataArray(
---> 74 ucat.convert_columns(
75 values, input_unit, 'kg m-2', p=psurf,
76 molar_mass=molar_mass
77 ),
78 dims=values.dims, attrs=values.attrs
79 )
81 if 'noise_level' in attrs:
82 noise_level = attrs['noise_level']
TypeError: convert_columns() got an unexpected keyword argument 'p'
[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)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, name)
1231 variable = self._variables[name]
1232 except KeyError:
-> 1233 _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
1234
KeyError: 'NO2_minus_estimated_background_mass'
During handling of the above exception, another exception occurred:
KeyError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py:1338, in Dataset.__getitem__(self, key)
1337 try:
-> 1338 return self._construct_dataarray(key)
1339 except KeyError as e:
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py:1233, in Dataset._construct_dataarray(self, name)
1232 except KeyError:
-> 1233 _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
1235 needed_dims = set(variable.dims)
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset_utils.py:79, in _get_virtual_variable(variables, key, dim_sizes)
78 if len(split_key) != 2:
---> 79 raise KeyError(key)
81 ref_name, var_name = split_key
KeyError: 'NO2_minus_estimated_background_mass'
The above exception was the direct cause of the following exception:
KeyError Traceback (most recent call last)
/tmp/ipykernel_999/1262715532.py in ?()
----> 1 results_csf = ddeq.csf.estimate_emissions(data, winds, sources, curves, 'NO2',
2 f_model=2.24, crs=CRS)
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/csf.py in ?(data, winds, sources, curves, gases, t_max, method, variable, crs, pixel_size, f_model, use_wind_timeseries)
882 current_method = 'sub-areas'
883 else:
884 current_method = 'gauss'
885
--> 886 ld = compute_line_density(
887 this, gases, variable,
888 pixel_size=pixel_size,
889 method=current_method,
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/csf.py in ?(data, gases, variable, xa, xb, ya, yb, dy, method, pixel_size, share_mu, share_sigma, background, extra_width)
504 polygon['along'] = xr.DataArray(0.5 * (xa + xb))
505
506 for gas in gases:
507 polygon.update(
--> 508 extract_pixels(data, gas, variable.format(gas=gas),
509 xa, xb, ya, yb, dy=dy)
510 )
511
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/csf.py in ?(data, gas, variable, xa, xb, ya, yb, dy)
402 dims='pixels')
403 polygon['y'] = xr.DataArray(data.yp.values[area], name='y',
404 dims='pixels')
405
--> 406 c = data[variable].values[area]
407 c[~isfinite] = np.nan
408 p = data['detected_plume'].values[area]
409
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, key)
1347
1348 # If someone attempts `ds['foo' , 'bar']` instead of `ds[['foo', 'bar']]`
1349 if isinstance(key, tuple):
1350 message += f"\nHint: use a list to select multiple variables, for example `ds[{list(key)}]`"
-> 1351 raise KeyError(message) from e
1352
1353 if utils.iterable_of_hashable(key):
1354 return self._copy_listed(key)
KeyError: "No variable named 'NO2_minus_estimated_background_mass'. Did you mean one of ('NO2_minus_estimated_background', 'NO2_estimated_background')?"
[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,
);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 fig = ddeq.vis.show_detected_plumes(
----> 2 data, curves, 1e6 * data['NO2'], gas='NO2', ld=results_csf,
3 vmin=0,
4 vmax=300,
5 label="NO$_2$ columns [µmol m$^{-2}$]",
6 winds=winds,
7 domain=DOMAIN,
8 do_zoom=False, show_clouds=False, crs=CRS,
9 figwidth=4,
10 );
NameError: name 'results_csf' is not defined
[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)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 2
1 fig, ax = plt.subplots(1, figsize=(8,2))
----> 2 ddeq.vis.plot_along_plume(ax, 'NO2', results_csf['Matimba'])
3 ax.set_xlim(right=207)
4 ax.set_ylim(0,150)
NameError: name 'results_csf' is not defined
[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()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 5
1 fig, axes = plt.subplots(1,3, figsize=(8,2), sharey=True)
4 for i, ax in zip([1,9,19], axes):
----> 5 r = results_csf['Matimba'].isel(along=i)
6 ddeq.vis.plot_across_section(r, gases=['NO2'], method='gauss', ax=ax,
7 legend='simple')
9 ax.text(-36, 19.5, '%d-%d km' % (r.xa/1e3, r.xb/1e3), ha='left', va='top')
NameError: name 'results_csf' is not defined
[10]:
ddeq.vis.plot_csf_result(['NO2'], data, winds, results_csf, curves, 'Matimba',
vmins=[0], vmaxs=[200e-6], crs=CRS);
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 ddeq.vis.plot_csf_result(['NO2'], data, winds, results_csf, curves, 'Matimba',
2 vmins=[0], vmaxs=[200e-6], crs=CRS);
NameError: name 'results_csf' is not defined
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)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, name)
1231 variable = self._variables[name]
1232 except KeyError:
-> 1233 _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
1234
KeyError: 'NO2_minus_estimated_background_mass'
During handling of the above exception, another exception occurred:
KeyError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py:1338, in Dataset.__getitem__(self, key)
1337 try:
-> 1338 return self._construct_dataarray(key)
1339 except KeyError as e:
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py:1233, in Dataset._construct_dataarray(self, name)
1232 except KeyError:
-> 1233 _, name, variable = _get_virtual_variable(self._variables, name, self.sizes)
1235 needed_dims = set(variable.dims)
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset_utils.py:79, in _get_virtual_variable(variables, key, dim_sizes)
78 if len(split_key) != 2:
---> 79 raise KeyError(key)
81 ref_name, var_name = split_key
KeyError: 'NO2_minus_estimated_background_mass'
The above exception was the direct cause of the following exception:
KeyError Traceback (most recent call last)
/tmp/ipykernel_999/2477038638.py in ?()
2 'NO2': {'Q': 3.0, # kg/s
3 'tau': 4*60**2 # seconds
4 }
5 }}
----> 6 data, results_gauss = ddeq.gauss.estimate_emissions(data, winds, sources, curves,
7 ['NO2'], priors=priors,
8 fit_decay_times=True)
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/gauss.py in ?(data, winds, sources, curves, gases, priors, variable, fit_decay_times, skip_overlapping_plumes, verbose)
651
652 this = data.sel(source=overlapping_sources)
653 wind_speeds = list(winds['speed'].sel(source=overlapping_sources).values)
654
--> 655 results, fields, CR = gaussian_plume_estimates(this, overlapping_sources,
656 curves, priors, wind_speeds,
657 variable, trace_gases=gases,
658 fit_decay_times=fit_decay_times,
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/gauss.py in ?(this, overlapping_sources, curves, priors, wind_speed, variable, trace_gases, fit_decay_times, sources, results, verbose)
412 if verbose:
413 params.pretty_print()
414
415 # Fit the data
--> 416 observations = this[variable.format(gas=gas)].values
417 mask = this.detected_plume.any('source')
418 mask = skimage.morphology.dilation(mask, skimage.morphology.square(5))
419 observations = np.where( mask, observations, np.nan )
~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/xarray/core/dataset.py in ?(self, key)
1347
1348 # If someone attempts `ds['foo' , 'bar']` instead of `ds[['foo', 'bar']]`
1349 if isinstance(key, tuple):
1350 message += f"\nHint: use a list to select multiple variables, for example `ds[{list(key)}]`"
-> 1351 raise KeyError(message) from e
1352
1353 if utils.iterable_of_hashable(key):
1354 return self._copy_listed(key)
KeyError: "No variable named 'NO2_minus_estimated_background_mass'. Did you mean one of ('NO2_minus_estimated_background', 'NO2_estimated_background')?"
[12]:
results_gauss = ddeq.emissions.convert_NO2_to_NOx_emissions(results_gauss, f=2.38)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 results_gauss = ddeq.emissions.convert_NO2_to_NOx_emissions(results_gauss, f=2.38)
NameError: name 'results_gauss' is not defined
[13]:
fig = ddeq.vis.plot_gauss_result(data, results_gauss, ['Matimba'],
'NO2', curves, crs=CRS, vmin=0, vmax=10e-6)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 fig = ddeq.vis.plot_gauss_result(data, results_gauss, ['Matimba'],
2 'NO2', curves, crs=CRS, vmin=0, vmax=10e-6)
NameError: name 'results_gauss' is not defined
[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())
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[15], line 1
----> 1 fig = ddeq.vis.show_detected_plumes(
2 data, curves,
3 1e6 * data['NO2_plume_model'].sel(source='Matimba'),
4 gas='NO2',
5 vmin=0,
6 vmax=300,
7 label="NO$_2$ column [µmol m$^{-2}$]",
8 winds=winds,
9 domain=DOMAIN,
10 do_zoom=False, show_clouds=False, crs=CRS,
11 figwidth=4, add_polygons=False
12 );
14 xp, yp = results_gauss[f'NO2_curve'].sel(source='Matimba').T
15 x, y = ddeq.misc.transform_coords(xp, yp, CRS, ccrs.PlateCarree(), use_xarray=False)
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/vis.py:905, in show_detected_plumes(data, curves, values, ld, add_polygons, zoom_on, figwidth, ax, cax, cmap, alpha, min_height, add_multiple_sources, vmin, vmax, show_clouds, winds, domain, add_origin, marker, markersize, label, do_zoom, crs, gas, dmin, delta, legend_loc, add_upstream_box, polygon_color, do_middle_lines)
902 curve = curves[source]
904 try:
--> 905 tmax = curve.arc2parameter(xb[-1])
906 except IndexError:
907 tmax = curve.t_o + 10e3
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/plume_coords.py:255, in Poly2D.arc2parameter(self, arc)
253 if self.interpolator is None:
254 ta = np.arange(self.tmin - 50e3, 2 * self.tmax, 500.0)
--> 255 la = [arc_length_of_2d_curve(self, self.t_o, v) for v in ta]
256 self.interpolator = scipy.interpolate.interp1d(la, ta)
258 return self.interpolator(arc)
File ~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/plume_coords.py:456, in arc_length_of_2d_curve(curve, a, b)
452 xt, yt = curve(t=t, m=1)
454 v = np.sqrt(xt**2 + yt**2)
--> 456 l = scipy.integrate.simps(v,t)
457 return l
AttributeError: module 'scipy.integrate' has no attribute 'simps'
