Example: Annual emissions#
Example of annual emissions using DIV and CSF with annual cycle for Berlin and Jänschwalde
[1]:
%matplotlib inline
import os
import itertools
import cartopy.crs as ccrs
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import ucat
import xarray as xr
import ddeq
CRS = ccrs.UTM(32)
WGS84 = ccrs.PlateCarree()
sources = ddeq.misc.read_point_sources()
CSF method#
The code below reads estimated CO\(_2\) and NO\(_X\) emissions from synthetic CO\(_2\) and NO\(_2\) observations generated from the SMARTCARB dataset. The results are provided in the data folder (SMARTCARB_CSF_example.nc).
The individual estimates are plotted for the Jänschwalde power plant and the city of Berlin. A seasonal cycle is fitted to the indivudal estimates and annual emissions are computed by temporal integration.
[2]:
d = xr.open_dataset(os.path.join(ddeq.DATA_PATH, 'SMARTCARB_CSF_example.nc'))
[3]:
fig, axes = plt.subplots(4,1, figsize=(7,7), sharex=True)
for i, (source, gas) in enumerate(itertools.product(['Janschwalde', 'Berlin'], ['CO2', 'NOx'])):
time = d['time']
units = 'Mt a$^{-1}$' if gas == 'CO2' else 'kt a$^{-1}$'
factor = 1e3 if gas == 'NOx' else 1.0
values = factor * d[f'{gas}_estimates'].sel(source=source)
values_std = factor * d[f'{gas}_estimates_std'].sel(source=source)
true = ddeq.smartcarb.read_true_emissions(None, 'NO2' if gas =='NOx' else gas, source)
true = true[true.index.hour == 11]
mask = (values > 0) & (values < 100) & (values_std < 20)
fit_result, func, times, cycle, chi2 = ddeq.timeseries.fit(
time.values[mask], values[mask], values_std[mask])
em, em_std = func.integrate(fit_result['x'], fit_result['x_std'])
axes[i].errorbar(time[mask].values, values[mask], yerr=values_std[mask],
ls='', marker='o', ms=4, label='Individual estimates')
axes[i].plot(true.index, true, '-', label=f'True emissions\n({np.mean(true):.1f} {units})')
axes[i].plot(times, cycle, 'k--', label=f'Estimated emissions\n({em:.1f}$\\pm${em_std:.1f} {units})')
axes[i].grid(True)
em, em_std = func.integrate(fit_result.x, fit_result.x_std)
axes[i].set_ylim(0,100)
axes[i].set_title(f'({"abcd"[i]}) {ddeq.vis.sub_numbers(gas)} emissions of {source}')
axes[i].set_ylabel(f'{ddeq.vis.sub_numbers(gas)} [{units}]')
axes[i].set_xlim(pd.Timestamp(2015,1,1), pd.Timestamp(2015,12,31))
axes[i].legend(ncol=1, loc=(1.02, 0.06))
axes[i].xaxis.set_major_formatter(mdates.DateFormatter('%b'))
axes[-1].set_xlabel('2015')
plt.tight_layout()
plt.savefig('CSF-timeseries.png', dpi=500)

Divergence method#
In this example, the divergence method was applied to one year of synthetic CO2M CO2 and NO2 images from the SMARTCARB dataset. The results are provided in the data folder (SMARTCARB_DIV_example.nc).
The code below can be used for reprocessing the example. Note that it requires the SMARTCARB dataset. The processing time is about 20 minutes.
[4]:
# on ICOS-CP computation time is about 1 minute per month and species
# (i.e. processing one year of CO2 and NO2 took about 20 minutes)
reprocess_div = False
if reprocess_div:
# change this to point to the SMARTCARB dataset for divergence method
ROOT = '/project/coco2/fileshare/WP4/SMARTCARB/'
SMARTCARB_DATA_PATH = os.path.join(ROOT, 'level2')
SMARTCARB_WIND_PATH = os.path.join(ROOT, 'wind_fields')
sources = ddeq.misc.read_point_sources()
sources = sources.sel(source=['Berlin', 'Janschwalde'])
l2 = ddeq.smartcarb.Level2Dataset(SMARTCARB_DATA_PATH, 'ace', co2_noise_scenario='low', no2_noise_scenario='low',
make_no2_error_cloud_dependent=False)
res = ddeq.div.estimate_emissions(
l2, SMARTCARB_WIND_PATH, sources=sources,
pattern='SMARTCARB_winds_%Y%m%d%H.nc',
wind_product='SMARTCARB', trace_gases=['CO2', 'NO2'],
varnames=['CO2', 'NO2'], smooth_data=[True, False],
remove_background=[True, False],
start_date='2015-01-01', end_date='2015-12-24'
)
res.to_netcdf('SMARTCARB_DIV_example.nc')
Open the example file and print estimated CO2 and NOX emissions:
[5]:
res = xr.open_dataset(os.path.join(ddeq.DATA_PATH, 'SMARTCARB_DIV_example.nc'))
annual = {
'CO2': {},
'NO2': {},
}
for source in ['Janschwalde', 'Berlin']:
for gas in ['CO2', 'NO2']:
em = res[f'{gas}_estimated_emissions'].sel(source=source)
em_std = res[f'{gas}_estimated_emissions_precision'].sel(source=source)
factor = 1.0 if gas == 'CO2' else 1.32e3
units = 'Mt' if gas == 'CO2' else 'kt'
em = factor * ddeq.misc.kgs_to_Mtyr(em)
em_std = factor * ddeq.misc.kgs_to_Mtyr(em_std)
print(f'{source:20s}{em:5.1f} ± {em_std:5.1f} {units} {gas} per year')
annual[gas][source] = f'{em:.1f} ± {em_std:.1f} {units}'
Janschwalde 43.2 ± 2.8 Mt CO2 per year
Janschwalde 37.3 ± 0.8 kt NO2 per year
Berlin 31.5 ± 2.8 Mt CO2 per year
Berlin 25.6 ± 0.1 kt NO2 per year
The following code plots the enhancements and fluxes of CO2 and NOX over Berlin and Jänschwalde from one year of SMARTCARB data:
[6]:
DOMAIN = ddeq.misc.Domain('BER+JAE', 12.5, 51.25, 15.3, 53.1)
[7]:
fig, axes = plt.subplots(2,2, figsize=(12,6.6),
subplot_kw={'projection': DOMAIN.proj})
# panel a
fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[0,0], admin_level=1, edgecolor='k')
for source in ['Berlin', 'Janschwalde']:
this = res.sel(source=source)
m = ax.pcolormesh(
this.longrid, this.latgrid,
ucat.convert_columns(this['CO2_grid'], 'm-2', 'ppmv', molar_mass='CO2'),
vmin=-0.4, vmax=0.4, transform=ccrs.PlateCarree())
plt.colorbar(m, ax=ax).set_label('XCO$_2$ enhancement [ppmv]')
# panel b
fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[0,1], admin_level=1, edgecolor='k')
for source in ['Berlin', 'Janschwalde']:
this = res.sel(source=source)
m = ax.pcolormesh(this.longrid, this.latgrid,
ucat.convert_columns(this['NO2_grid'], 'm-2', 'umol m-2', molar_mass='NO2'),
vmin=0, vmax=100, transform=ccrs.PlateCarree())
plt.colorbar(m, ax=ax).set_label('NO$_2$ enhancement [µmol m$^{-2}$]')
# panel c
fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[1,0], admin_level=1, edgecolor='k')
for source in ['Berlin', 'Janschwalde']:
this = res.sel(source=source)
m = ax.pcolormesh(this.longrid, this.latgrid,
this['CO2_div'] / 1e3,
vmin=-2, vmax=4, transform=ccrs.PlateCarree())
plt.colorbar(m, ax=ax).set_label('XCO$_2$ flux [kg m$^{-2}$ s$^{-1}$]')
ax.text(0.975, 0.96, f'Berlin: {annual["CO2"]["Berlin"]}\nJänschwalde: {annual["CO2"]["Janschwalde"]}',
ha='right', va='top', fontsize='small', transform=ax.transAxes,
bbox={'ec': 'black', 'fc': 'white', 'lw': 0.15}
)
# panel d
fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[1,1], admin_level=1, edgecolor='k')
for source in ['Berlin', 'Janschwalde']:
this = res.sel(source=source)
m = ax.pcolormesh(this.longrid, this.latgrid,
this['NO2_div'],
vmin=-2, vmax=4, transform=ccrs.PlateCarree())
plt.colorbar(m, ax=ax).set_label('NO$_x$ flux [g m$^{-2}$ s$^{-1}$]')
ax.text(0.975, 0.96, f'Berlin: {annual["NO2"]["Berlin"]}\nJänschwalde: {annual["NO2"]["Janschwalde"]}',
ha='right', va='top', fontsize='small', transform=ax.transAxes,
bbox={'ec': 'black', 'fc': 'white', 'lw': 0.15}
)
for i, ax in enumerate(axes.flatten()):
ddeq.vis.add_gridlines(ax, dlon=0.5, dlat=0.5)
ddeq.vis.add_hot_spots(
ax,
sources=sources.sel(source=['Berlin', 'Janschwalde', 'Schwarze Pumpe', 'Boxberg']),
do_path_effect=True, ms=3, size='small'
)
ax.set_title(f'({"abcd"[i]})')
plt.tight_layout()
plt.savefig('DIV-example.png', dpi=500)
