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
[2]:
<xarray.Dataset>
Dimensions:  (source: 22)
Coordinates:
  * source   (source) object 'Belchatow' 'Berlin' ... 'Staudinger' 'Turow'
Data variables:
    lon_o    (source) float64 19.33 13.41 14.57 15.44 ... 14.36 8.957 14.91
    lat_o    (source) float64 51.27 52.52 51.42 50.02 ... 51.53 50.09 50.94
    type     (source) object 'power plant' 'city' ... 'power plant'
    label    (source) object 'Bełchatów' 'Berlin' ... 'Staudinger' 'Turów'
    radius   (source) float64 5.0 15.0 5.0 5.0 5.0 15.0 ... 5.0 5.0 5.0 5.0 5.0

Synthetic satellite observations#

Synthetic satellite observations are available from the SMARTCARB project (https://doi.org/10.5281/zenodo.4048227). The ddeq package can read the data files and automatically applies random noise and cloud filters to the observations. The code also fixes some issues with the dataset such as wrong emissions for industry in Berlin in January and July. It is also possible to scale the anthropogenic model tracers:

[3]:
filename = os.path.join(DATA_PATH,'Sentinel_7_CO2_2015042311_o1670_l0483.nc')
data_level2 = ddeq.smartcarb.read_level2(filename, co2_noise_scenario='medium',
                                         no2_noise_scenario='high')
data_level2
[3]:
<xarray.Dataset>
Dimensions:  (nobs: 811, nrows: 123, ncorners: 4)
Dimensions without coordinates: nobs, nrows, ncorners
Data variables:
    time     datetime64[ns] 2015-04-23T11:00:00
    lon      (nobs, nrows) float32 16.18 16.21 16.25 16.29 ... 13.3 13.32 13.35
    lat      (nobs, nrows) float32 60.54 60.54 60.53 60.53 ... 46.03 46.03 46.02
    lonc     (nobs, nrows, ncorners) float32 16.16 16.2 16.19 ... 13.36 13.33
    latc     (nobs, nrows, ncorners) float32 60.55 60.55 60.53 ... 46.01 46.02
    clouds   (nobs, nrows) float32 nan nan nan nan nan ... nan nan nan nan nan
    psurf    (nobs, nrows) float32 nan nan nan nan nan ... nan nan nan nan nan
    CO2      (nobs, nrows) float32 nan nan nan nan nan ... nan nan nan nan nan
    CO2_std  (nobs, nrows) float32 0.7 0.7 0.7 0.7 0.7 ... 0.7 0.7 0.7 0.7 0.7
    NO2      (nobs, nrows) float32 nan nan nan nan nan ... nan nan nan nan nan
    NO2_std  (nobs, nrows) float32 2e+15 2e+15 2e+15 2e+15 ... 2e+15 2e+15 2e+15
Attributes:
    satellite:    CO2M
    orbit:        1670
    lon_eq:       483
    DESCRIPTION:  Synthetic XCO2 and NO2 satellite image with auxiliary data ...
    DATAORIGIN:   SMARTCARB study
    DOI:          10.5281/zenodo.4048227
    CREATOR:      Gerrit Kuhlmann et al.
    EMAIL:        gerrit.kuhlmann@empa.ch
    AFFILIATION:  Empa Duebendorf, Switzerland

The data can easily be plotted using the ddeq.vis.visualize function, which requires the satellite data, the name of the trace gas, the SMARTCARB model domain, and the dataset of sources for labeling point sources.

[4]:
fig = ddeq.vis.visualize(data_level2, data_level2.NO2, gas='NO2', domain=DOMAIN,
                         sources=sources)
_images/tutorial-introduction-to-ddeq_8_0.png

It is also possible to read and visualize the (vertically integrated) COSMO-GHG fields:

[5]:
time = pd.Timestamp(data_level2.time.to_pandas())
filename = os.path.join(DATA_PATH, 'cosmo_2d_2015042311.nc')

data_cosmo = ddeq.smartcarb.read_cosmo(filename, 'CO2')

fig = ddeq.vis.make_field_map(data_cosmo, trace_gas='CO2', domain=DOMAIN,
                              vmin=404, vmax=408, alpha=data_cosmo['CLCT'],
                              border=50.0, label='XCO$_2$ [ppm]')
_images/tutorial-introduction-to-ddeq_10_0.png

Exercise#

  • Read SMARTCARB Level-2 data for 23 April 2015, 11 UTC (orbit: 1670, lon_eq: 0483) using a low-noise CO2 and high-noise NO2 uncertainty scenario.

  • Plot the XCO2 observations using ddeq.vis.visualize, mask cloud fractions larger than 1%, and add labels point sources (Berlin, Boxberg, Janschwalde, Lippendorf, Schwarze Pumpe and Turow).

  • Read and add the XCO2 field from the COSMO-GHG model (ddeq.smartcarb.read_trace_gas_field) and additional fields (ddeq.smartcarb.read_fields) to the plot.

  • Add a square showing the study area given by lower left and upper right points of 12.0°N, 50.7°E and 15.5°N, 52.7°N, respectively.

[ ]:

A solution can be found in the example-hakkarainen-2022-fig-01.ipynb file.

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