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 matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import pandas as pd
import xarray as xr

# import and setup matplotlib
%matplotlib inline
plt.rcParams['figure.dpi'] = 100

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

CRS = ccrs.UTM(32)
WGS84 = ccrs.PlateCarree()

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

[2]:

# list of point sources
sources = ddeq.misc.read_point_sources()
sources

Wind fields#

Data-driven emission quantification always requires a wind speed to convert (integrated) enhancements to fluxes. ddeq.wind provides access to different wind datasets such as ERA-5 and the SMARTCARB dataset. The example below returns the mean winds within a 0.05° radius around each source in sources from the SMARTCARB dataset:

[6]:

winds = ddeq.wind.read_at_sources(time, sources, radius=0.05,
                                  product='SMARTCARB', data_path=DATA_PATH)

It is also possible to show the wind direction at the source location providing winds as an argument to ddeq.vis.visualize:

[7]:

fig = ddeq.vis.visualize(data_level2, data_level2.NO2, gas='NO2', domain=DOMAIN,
                         sources=sources, winds=winds)

_images/tutorial-introduction-to-ddeq_18_0.png

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)

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.visualize function:

[10]:

fig = ddeq.vis.visualize(data, data['NO2'], gas='NO2', domain=DOMAIN,
                         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.visualize

[ ]:

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]:

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
)

data = ddeq.dplume.detect_plumes(
    data, sources.sel(source=['Berlin', 'Janschwalde']),
    variable='NO2', variable_std='NO2_std',
    filter_type='gaussian', filter_size=0.5
)

# setting `do_zoom` to True will zoom on detected plumes and `draw_gridlines`
# will add lines
ddeq.vis.visualize(
    data, 'NO2', gas='NO2', domain=DOMAIN, winds=winds, do_zoom=True,
    show_clouds=True, draw_gridlines=True
);

_images/tutorial-introduction-to-ddeq_29_0.png

[12]:

data, curves = ddeq.plume_coords.compute_plume_line_and_coords(
    data, crs=CRS, radius=25e3, plume_area='area'
)

The following code shows the result:

[13]:

ddeq.vis.show_detected_plumes(
    data, curves, data['NO2'], gas='NO2', domain=DOMAIN, winds=winds,
    do_zoom=True, crs=CRS
);

_images/tutorial-introduction-to-ddeq_32_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.

[14]:

data = ddeq.emissions.prepare_data(data, 'CO2')
data = ddeq.emissions.prepare_data(data, 'NO2')

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

[15]:

ddeq.vis.visualize(
    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}$]'
);

_images/tutorial-introduction-to-ddeq_37_0.png

Cross-sectional flux method#

The following code estimated CO\(_2\) and NO\(_x\) emissions for a point source (Jänschwalde) and a city (Berlin):

Cross-sectional flux method for point source#

Wind speed and direction are taken from ERA-5 data that are downloaded from the Copernicus Climate Data Store (CDS). This can be done automatically using cdsapi but requires a CDS account and it might be slow especially when downloading ERA-5 on model levels.

In the example below, winds are computed from ERA-5 on model levels using the GNFR-A emission profile for vertical averaging. A subset of ERA-5 data from the SMARTCARB model domain is included in DATA_PATH for testing.

[16]:

winds = ddeq.wind.read_at_sources(data.time, data, product='ERA5', era5_prefix='SMARTCARB',
                                  data_path=DATA_PATH)

The cross-sectional flux (csf) method is performed by the following function.

Note that f_model gives the factor for converting NO2 to NOx line densities.

[17]:

results = ddeq.csf.estimate_emissions(data, winds, sources, curves,
                                      method='gauss', gases=['CO2', 'NO2'], crs=CRS,
                                      f_model=1.32
                                     )

The results can be visualized with the following function:

[18]:

with xr.set_options(keep_attrs=True):
    fig = ddeq.vis.plot_csf_result(
        ['CO2', 'NO2'],
        data, winds, results,
        source='Janschwalde',
        curves=curves,
        domain=DOMAIN, crs=CRS
    )

plt.savefig('ddeq-example.png',  bbox_inches='tight')

_images/tutorial-introduction-to-ddeq_45_0.png

Cross-sectional flux method for a city#

The cross sectional flux method over a city is slightly different, because the flux slowly builds up over the city area. For NO\(_x\), fluxes over the city are therefore modeled by a Gaussian curve.

[19]:

fig = ddeq.vis.plot_csf_result(
    ['CO2', 'NO2'],
    data, winds, results,
    source='Berlin',
    curves=curves,
    domain=DOMAIN, crs=CRS,
)

_images/tutorial-introduction-to-ddeq_47_0.png

Exercise#

Write code to quantify the CO2 and NOx emissions of other point sources.

[ ]:

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_at_sources(data.time, data, product='SMARTCARB',
                                  data_path=DATA_PATH)

[21]:

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

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

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

                Estimate        True
Berlin          17.7 Mt/a       23.3 Mt/a
Janschwalde     38.5 Mt/a       42.4 Mt/a

[22]:

results = ddeq.ime.estimate_emissions(data, winds, sources, curves,
                                      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_estimated_emissions'].sel(source=name).values
        Q_true = ddeq.smartcarb.read_true_emissions(
            pd.Timestamp(data.time.values), 'NO2', name
        ).mean()
        print(f'{name:10s}\t{ddeq.misc.kgs_to_Mtyr(Q*1000):.1f} kt/a\t{Q_true:.1f} kt/a')

                Estimate        True
Berlin          32.7 kt/a       24.7 kt/a
Janschwalde     42.0 kt/a       34.2 kt/a

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#

Cross-sectional flux method for point source#

Cross-sectional flux method for a city#

Exercise#

Integrated mass enhancement#