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'
)
_images/example_tropomi_Matimba_3_0.png

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,
);
_images/example_tropomi_Matimba_8_0.png
[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>
_images/example_tropomi_Matimba_9_1.png
[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()

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

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)
_images/example_tropomi_Matimba_15_0.png
[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>]
_images/example_tropomi_Matimba_17_1.png