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

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.

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

<xarray.Dataset> Size: 1kB
Dimensions:  (source: 22)
Coordinates:
  * source   (source) object 176B 'Belchatow' 'Berlin' ... 'Staudinger' 'Turow'
Data variables:
    lon_o    (source) float64 176B 19.33 13.41 14.57 15.44 ... 14.36 8.957 14.91
    lat_o    (source) float64 176B 51.27 52.52 51.42 50.02 ... 51.53 50.09 50.94
    type     (source) object 176B 'power plant' 'city' ... 'power plant'
    label    (source) object 176B 'Bełchatów' 'Berlin' ... 'Staudinger' 'Turów'
    radius   (source) float64 176B 5.0 15.0 5.0 5.0 5.0 ... 5.0 5.0 5.0 5.0 5.0

xarray.Dataset

• Dimensions:
  • source: 22
• Coordinates: (1)
  • source
    (source)
    object
    'Belchatow' 'Berlin' ... 'Turow'

    array(['Belchatow', 'Berlin', 'Boxberg', 'Chvaletice', 'Dolna Odra', 'Hamburg',
           'Heyden', 'Janschwalde', 'Lipetsk', 'Lippendorf', 'Matimba', 'Melnik',
           'Opole', 'Patnow', 'Paris', 'Pocerady', 'Prague', 'Prunerov',
           'Schkopau', 'Schwarze Pumpe', 'Staudinger', 'Turow'], dtype=object)

• Data variables: (5)
  • lon_o
    (source)
    float64
    19.33 13.41 14.57 ... 8.957 14.91

    name :

    longitude of point source

    array([19.3267    , 13.407696  , 14.57469559, 15.4437438 , 14.466333  ,
           10.001389  ,  8.98862953, 14.45349026, 39.6196    , 12.37124538,
           27.610556  , 14.41942402, 17.88227332, 18.23756027,  2.352222  ,
           13.68334331, 14.416667  , 13.25937817, 11.94972392, 14.36109257,
            8.95729041, 14.91128159])

  • lat_o
    (source)
    float64
    51.27 52.52 51.42 ... 50.09 50.94

    name :

    latitude of point source

    array([ 51.2659    ,  52.516984  ,  51.41649246,  50.01838901,
            53.205776  ,  53.565278  ,  52.38599786,  51.84154511,
            52.5603    ,  51.18745041, -23.668333  ,  50.41045529,
            50.74892801,  52.30444336,  48.856667  ,  50.42475412,
            50.083333  ,  50.41673025,  51.39481843,  51.5342865 ,
            50.08554139,  50.94282532])

  • type
    (source)
    object
    'power plant' ... 'power plant'

    name :

    type of point source

    array(['power plant', 'city', 'power plant', 'power plant', 'power plant',
           'city', 'power plant', 'power plant', 'steel plant', 'power plant',
           'power plant', 'power plant', 'power plant', 'power plant', 'city',
           'power plant', 'city', 'power plant', 'power plant', 'power plant',
           'power plant', 'power plant'], dtype=object)

  • label
    (source)
    object
    'Bełchatów' 'Berlin' ... 'Turów'

    array(['Bełchatów', 'Berlin', 'Boxberg', 'Chvaletice', 'Dolna Odra',
           'Hamburg', 'Heyden', 'Jänschwalde', 'Lipetsk', 'Lippendorf',
           'Matimba & Medupi', 'Mělník', 'Opole', 'Patnów', 'Paris',
           'Počerady', 'Prague', 'Prunéřov', 'Schkopau', 'Schwarze Pumpe',
           'Staudinger', 'Turów'], dtype=object)

  • radius
    (source)
    float64
    5.0 15.0 5.0 5.0 ... 5.0 5.0 5.0

    name :

    Radius of point source around location

    units :

    km

    array([ 5., 15.,  5.,  5.,  5., 15.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,
            5., 15.,  5., 15.,  5.,  5.,  5.,  5.,  5.])

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:

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

<xarray.Dataset> Size: 7MB
Dimensions:  (nobs: 811, nrows: 123, ncorners: 4)
Dimensions without coordinates: nobs, nrows, ncorners
Data variables:
    time     datetime64[ns] 8B 2015-04-23T11:00:00
    lon      (nobs, nrows) float32 399kB 16.18 16.21 16.25 ... 13.3 13.32 13.35
    lat      (nobs, nrows) float32 399kB 60.54 60.54 60.53 ... 46.03 46.03 46.02
    lonc     (nobs, nrows, ncorners) float32 2MB 16.16 16.2 ... 13.36 13.33
    latc     (nobs, nrows, ncorners) float32 2MB 60.55 60.55 ... 46.01 46.02
    clouds   (nobs, nrows) float32 399kB nan nan nan nan nan ... nan nan nan nan
    psurf    (nobs, nrows) float32 399kB nan nan nan nan nan ... nan nan nan nan
    CO2      (nobs, nrows) float32 399kB nan nan nan nan nan ... nan nan nan nan
    CO2_std  (nobs, nrows) float32 399kB 0.7 0.7 0.7 0.7 0.7 ... 0.7 0.7 0.7 0.7
    NO2      (nobs, nrows) float64 798kB nan nan nan nan nan ... nan nan nan nan
    NO2_std  (nobs, nrows) float32 399kB 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

xarray.Dataset

• Dimensions:
  • nobs: 811
  • nrows: 123
  • ncorners: 4
• Data variables: (11)
  • time
    ()
    datetime64[ns]
    2015-04-23T11:00:00

    array('2015-04-23T11:00:00.000000000', dtype='datetime64[ns]')

  • lon
    (nobs, nrows)
    float32
    16.18 16.21 16.25 ... 13.32 13.35

    standard_name :

    longitude

    long_name :

    pixel center longitude

    units :

    degrees east

    array([[16.177404, 16.214611, 16.251749, ..., 20.429628, 20.465448,
            20.501312],
           [16.16753 , 16.204718, 16.241838, ..., 20.417713, 20.453516,
            20.489363],
           [16.157665, 16.194834, 16.231936, ..., 20.405807, 20.441595,
            20.477425],
           ...,
           [10.197353, 10.224309, 10.25122 , ..., 13.308682, 13.335179,
            13.361712],
           [10.191625, 10.218572, 10.245474, ..., 13.302004, 13.328492,
            13.355018],
           [10.185898, 10.212837, 10.239731, ..., 13.295328, 13.321809,
            13.348327]], shape=(811, 123), dtype=float32)

  • lat
    (nobs, nrows)
    float32
    60.54 60.54 60.53 ... 46.03 46.02

    standard_name :

    latitude

    long_name :

    pixel center latitude

    units :

    degrees north

    array([[60.543377, 60.538837, 60.534294, ..., 59.947453, 59.941753,
            59.936035],
           [60.526134, 60.521595, 60.517056, ..., 59.93052 , 59.924824,
            59.91911 ],
           [60.508892, 60.504356, 60.499817, ..., 59.913586, 59.907894,
            59.902184],
           ...,
           [46.493587, 46.490265, 46.486946, ..., 46.06552 , 46.06148 ,
            46.05743 ],
           [46.4761  , 46.472782, 46.46946 , ..., 46.04817 , 46.044132,
            46.04008 ],
           [46.458614, 46.455296, 46.451977, ..., 46.03082 , 46.026783,
            46.022736]], shape=(811, 123), dtype=float32)

  • lonc
    (nobs, nrows, ncorners)
    float32
    16.16 16.2 16.19 ... 13.36 13.33

    standard_name :

    longitude

    long_name :

    pixel corners longitude

    units :

    degrees east

    array([[[16.163712 , 16.200962 , 16.191074 , 16.153843 ],
            [16.200962 , 16.238144 , 16.228235 , 16.191074 ],
            [16.238144 , 16.275257 , 16.265331 , 16.228235 ],
            ...,
            [20.417692 , 20.4535   , 20.44157  , 20.405779 ],
            [20.4535   , 20.48935  , 20.477404 , 20.44157  ],
            [20.48935  , 20.525242 , 20.513279 , 20.477404 ]],
    
           [[16.153843 , 16.191074 , 16.181194 , 16.143982 ],
            [16.191074 , 16.228235 , 16.218338 , 16.181194 ],
            [16.228235 , 16.265331 , 16.255415 , 16.218338 ],
            ...,
            [20.405779 , 20.44157  , 20.429651 , 20.393877 ],
            [20.44157  , 20.477404 , 20.465467 , 20.429651 ],
            [20.477404 , 20.513279 , 20.501328 , 20.465467 ]],
    
           [[16.143982 , 16.181194 , 16.171324 , 16.13413  ],
            [16.181194 , 16.218338 , 16.20845  , 16.171324 ],
            [16.218338 , 16.255415 , 16.245508 , 16.20845  ],
            ...,
    ...
            [13.298786 , 13.325269 , 13.318584 , 13.2921095],
            [13.325269 , 13.351788 , 13.345096 , 13.318584 ],
            [13.351788 , 13.378345 , 13.371644 , 13.345096 ]],
    
           [[10.180996 , 10.20797  , 10.202238 , 10.175272 ],
            [10.20797  , 10.2348995, 10.229158 , 10.202238 ],
            [10.2348995, 10.261784 , 10.256035 , 10.229158 ],
            ...,
            [13.2921095, 13.318584 , 13.311903 , 13.285437 ],
            [13.318584 , 13.345096 , 13.338407 , 13.311903 ],
            [13.345096 , 13.371644 , 13.364947 , 13.338407 ]],
    
           [[10.175272 , 10.202238 , 10.196509 , 10.169551 ],
            [10.202238 , 10.229158 , 10.223421 , 10.196509 ],
            [10.229158 , 10.256035 , 10.250289 , 10.223421 ],
            ...,
            [13.285437 , 13.311903 , 13.305225 , 13.278768 ],
            [13.311903 , 13.338407 , 13.331721 , 13.305225 ],
            [13.338407 , 13.364947 , 13.358254 , 13.331721 ]]],
          shape=(811, 123, 4), dtype=float32)

  • latc
    (nobs, nrows, ncorners)
    float32
    60.55 60.55 60.53 ... 46.01 46.02

    standard_name :

    latitude

    long_name :

    pixel corners latitude

    units :

    degrees north

    array([[[60.554264, 60.54973 , 60.532486, 60.53702 ],
            [60.54973 , 60.545185, 60.527946, 60.532486],
            [60.545185, 60.540638, 60.523403, 60.527946],
            ...,
            [59.95876 , 59.953068, 59.93614 , 59.941826],
            [59.953068, 59.947357, 59.93043 , 59.93614 ],
            [59.947357, 59.94163 , 59.924706, 59.93043 ]],
    
           [[60.53702 , 60.532486, 60.515244, 60.51978 ],
            [60.532486, 60.527946, 60.510708, 60.515244],
            [60.527946, 60.523403, 60.50617 , 60.510708],
            ...,
            [59.941826, 59.93614 , 59.91921 , 59.924892],
            [59.93614 , 59.93043 , 59.913506, 59.91921 ],
            [59.93043 , 59.924706, 59.907784, 59.913506]],
    
           [[60.51978 , 60.515244, 60.498   , 60.502533],
            [60.515244, 60.510708, 60.49347 , 60.498   ],
            [60.510708, 60.50617 , 60.48893 , 60.49347 ],
            ...,
    ...
            [46.07621 , 46.072178, 46.05483 , 46.05886 ],
            [46.072178, 46.06813 , 46.050785, 46.05483 ],
            [46.06813 , 46.064075, 46.046726, 46.050785]],
    
           [[46.486504, 46.483185, 46.4657  , 46.469013],
            [46.483185, 46.479862, 46.46238 , 46.4657  ],
            [46.479862, 46.476543, 46.45906 , 46.46238 ],
            ...,
            [46.05886 , 46.05483 , 46.03748 , 46.041508],
            [46.05483 , 46.050785, 46.033436, 46.03748 ],
            [46.050785, 46.046726, 46.02938 , 46.033436]],
    
           [[46.469013, 46.4657  , 46.44821 , 46.451527],
            [46.4657  , 46.46238 , 46.444893, 46.44821 ],
            [46.46238 , 46.45906 , 46.441574, 46.444893],
            ...,
            [46.041508, 46.03748 , 46.020126, 46.024155],
            [46.03748 , 46.033436, 46.016087, 46.020126],
            [46.033436, 46.02938 , 46.01203 , 46.016087]]],
          shape=(811, 123, 4), dtype=float32)

  • clouds
    (nobs, nrows)
    float32
    nan nan nan nan ... nan nan nan nan

    standard_name :

    cloud_area_fraction

    units :

    1

    long_name :

    total cloud cover

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]],
          shape=(811, 123), dtype=float32)

  • psurf
    (nobs, nrows)
    float32
    nan nan nan nan ... nan nan nan nan

    standard_name :

    surface_air_pressure

    units :

    Pa

    long_name :

    surface pressure

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]],
          shape=(811, 123), dtype=float32)

  • CO2
    (nobs, nrows)
    float32
    nan nan nan nan ... nan nan nan nan

    VEG50_uncertainties :

    0.5, 0.7, 1.0

    standard_name :

    XCO2 random uncertainty using error parametrization formula

    scaling_factors :

    0.333, 0.500, 0.667

    long_name :

    XCO2 random uncertainty using error parametrization formula

    noise_level :

    0.7

    cloud_threshold :

    0.01

    units :

    ppmv

    standard name :

    CO2_column-averaged_dry-air_mole_fraction with random noise

    long name :

    CO2_column-averaged_dry-air_mole_fraction with random noise

    scaling :

    1.0

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]],
          shape=(811, 123), dtype=float32)

  • CO2_std
    (nobs, nrows)
    float32
    0.7 0.7 0.7 0.7 ... 0.7 0.7 0.7 0.7

    units :

    ppmv

    standard name :

    CO2_column-averaged_dry-air_mole_fraction

    long name :

    CO2_column-averaged_dry-air_mole_fraction

    array([[0.7, 0.7, 0.7, ..., 0.7, 0.7, 0.7],
           [0.7, 0.7, 0.7, ..., 0.7, 0.7, 0.7],
           [0.7, 0.7, 0.7, ..., 0.7, 0.7, 0.7],
           ...,
           [0.7, 0.7, 0.7, ..., 0.7, 0.7, 0.7],
           [0.7, 0.7, 0.7, ..., 0.7, 0.7, 0.7],
           [0.7, 0.7, 0.7, ..., 0.7, 0.7, 0.7]],
          shape=(811, 123), dtype=float32)

  • NO2
    (nobs, nrows)
    float64
    nan nan nan nan ... nan nan nan nan

    noise_level :

    2000000000000000.0

    cloud_threshold :

    0.3

    units :

    cm-2

    standard name :

    NO2_vertical_column_density with random noise

    long name :

    NO2_vertical_column_density with random noise

    scaling :

    1.0

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], shape=(811, 123))

  • NO2_std
    (nobs, nrows)
    float32
    2e+15 2e+15 2e+15 ... 2e+15 2e+15

    units :

    cm-2

    standard name :

    NO2_vertical_column_density

    long name :

    NO2_vertical_column_density

    array([[2.e+15, 2.e+15, 2.e+15, ..., 2.e+15, 2.e+15, 2.e+15],
           [2.e+15, 2.e+15, 2.e+15, ..., 2.e+15, 2.e+15, 2.e+15],
           [2.e+15, 2.e+15, 2.e+15, ..., 2.e+15, 2.e+15, 2.e+15],
           ...,
           [2.e+15, 2.e+15, 2.e+15, ..., 2.e+15, 2.e+15, 2.e+15],
           [2.e+15, 2.e+15, 2.e+15, ..., 2.e+15, 2.e+15, 2.e+15],
           [2.e+15, 2.e+15, 2.e+15, ..., 2.e+15, 2.e+15, 2.e+15]],
          shape=(811, 123), dtype=float32)

• Attributes: (9)

  satellite :

  CO2M

  orbit :

  1670

  lon_eq :

  483

  DESCRIPTION :

  Synthetic XCO2 and NO2 satellite image with auxiliary data for estimating CO2/NOx emissions of cities and point sources

  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.

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

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

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]')

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:

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:

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

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).

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

fig = ddeq.vis.visualize(data, data['NO2'], gas='NO2', domain=DOMAIN,
                         winds=winds, do_zoom=False, show_clouds=True)

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.

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
);

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

The following code shows the result:

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

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[13], line 1
----> 1 ddeq.vis.show_detected_plumes(
      2     data, curves, data['NO2'], gas='NO2', domain=DOMAIN, winds=winds,
      3     do_zoom=True, crs=CRS
      4 );

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'

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.

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

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[14], line 1
----> 1 data = ddeq.emissions.prepare_data(data, 'CO2')
      2 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'

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

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}$]'
);

---------------------------------------------------------------------------
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: 'CO2_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: 'CO2_minus_estimated_background_mass'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_1052/237838958.py in ?()
----> 1 ddeq.vis.visualize(
      2     data, 'CO2_minus_estimated_background_mass', gas='CO2', domain=DOMAIN,
      3     winds=None, do_zoom=True, show_clouds=True, draw_gridlines=True,
      4     vmin=-20e-3, vmax=40e-3,

~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/vis.py in ?(data, values, do_zoom, winds, show_clouds, alpha, vlim, legend_loc, simple_legend, add_multiple_sources, gas, scaling, sources, names, names_align, figwidth, vmin, vmax, zoom_on, min_height, ax, cax, cmap, marker, markersize, domain, label, units, draw_gridlines)
    566     else:
    567         fig = ax, cax
    568
    569     if isinstance(values, str):
--> 570         values = data[values]
    571
    572     if scaling is not None:
    573         with xr.set_options(keep_attrs=True):

~/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 'CO2_minus_estimated_background_mass'. Did you mean one of ('CO2_minus_estimated_background', 'CO2_estimated_background')?"

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.

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.

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

---------------------------------------------------------------------------
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: 'CO2_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: 'CO2_minus_estimated_background_mass'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_1052/3884418323.py in ?()
----> 1 results = ddeq.csf.estimate_emissions(data, winds, sources, curves,
      2                                       method='gauss', gases=['CO2', 'NO2'], crs=CRS,
      3                                       f_model=1.32
      4                                      )

~/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 'CO2_minus_estimated_background_mass'. Did you mean one of ('CO2_minus_estimated_background', 'CO2_estimated_background')?"

The results can be visualized with the following function:

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')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[18], line 4
      1 with xr.set_options(keep_attrs=True):
      2     fig = ddeq.vis.plot_csf_result(
      3         ['CO2', 'NO2'],
----> 4         data, winds, results,
      5         source='Janschwalde',
      6         curves=curves,
      7         domain=DOMAIN, crs=CRS
      8     )
     10 plt.savefig('ddeq-example.png',  bbox_inches='tight')

NameError: name 'results' is not defined

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.

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[19], line 3
      1 fig = ddeq.vis.plot_csf_result(
      2     ['CO2', 'NO2'],
----> 3     data, winds, results,
      4     source='Berlin',
      5     curves=curves,
      6     domain=DOMAIN, crs=CRS,
      7 )

NameError: name 'results' is not defined

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.

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

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')

---------------------------------------------------------------------------
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: 'CO2_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: 'CO2_minus_estimated_background_mass'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_1052/2522278214.py in ?()
----> 1 results = ddeq.ime.estimate_emissions(data, winds, sources, curves, gas='CO2')
      2
      3 print(' ' * 10, '\tEstimate\tTrue')
      4

~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/ime.py in ?(data, winds, sources, curves, gas, variable, L_min, L_max, decay_time, min_pixel_number)
    100         # plume area
    101         area = this['plume_area'].values
    102
    103         # interpolate missing data
--> 104         vcd = this[variable.format(gas=gas)].values
    105         missing_data = np.isfinite(vcd)
    106
    107         # fill cloud- and data-gaps

~/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 'CO2_minus_estimated_background_mass'. Did you mean one of ('CO2_minus_estimated_background', 'CO2_estimated_background')?"

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')

---------------------------------------------------------------------------
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_1052/1792275635.py in ?()
----> 1 results = ddeq.ime.estimate_emissions(data, winds, sources, curves,
      2                                       gas='NO2', decay_time=4*60**2)
      3
      4 results = ddeq.emissions.convert_NO2_to_NOx_emissions(results, f=1.32)

~/checkouts/readthedocs.org/user_builds/ddeq/envs/v1.0/lib/python3.12/site-packages/ddeq/ime.py in ?(data, winds, sources, curves, gas, variable, L_min, L_max, decay_time, min_pixel_number)
    100         # plume area
    101         area = this['plume_area'].values
    102
    103         # interpolate missing data
--> 104         vcd = this[variable.format(gas=gas)].values
    105         missing_data = np.isfinite(vcd)
    106
    107         # fill cloud- and data-gaps

~/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 ('CO2_minus_estimated_background', 'CO2_estimated_background')?"