{ "cells": [ { "cell_type": "markdown", "id": "4cd53a91-b803-4750-bb29-64fe2bbcfda1", "metadata": {}, "source": [ "# Example: Annual emissions" ] }, { "cell_type": "markdown", "id": "11816043-7f96-4170-b039-aee4fcb5df17", "metadata": {}, "source": [ "Example of annual emissions using DIV and CSF with annual cycle for Berlin and Jänschwalde" ] }, { "cell_type": "code", "execution_count": null, "id": "190736d1-fa5c-4af0-a3e0-4446137d6e05", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import os\n", "import itertools\n", "\n", "import cartopy.crs as ccrs\n", "import matplotlib.dates as mdates\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import ucat\n", "import xarray as xr\n", "\n", "import ddeq\n", "\n", "CRS = ccrs.UTM(32)\n", "WGS84 = ccrs.PlateCarree()\n", "\n", "sources = ddeq.misc.read_point_sources()" ] }, { "cell_type": "markdown", "id": "7d53ad73-6b73-4fd6-a29c-839f43fb128a", "metadata": {}, "source": [ "## CSF method\n", "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).\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "34ee4ee3-0f93-4491-93f0-5b66f0a68927", "metadata": {}, "outputs": [], "source": [ "d = xr.open_dataset(os.path.join(ddeq.DATA_PATH, 'SMARTCARB_CSF_example.nc'))" ] }, { "cell_type": "code", "execution_count": null, "id": "5a5074dc-9a59-4ab8-af60-61eeaac0b5dc", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(4,1, figsize=(7,7), sharex=True)\n", "\n", "for i, (source, gas) in enumerate(itertools.product(['Janschwalde', 'Berlin'], ['CO2', 'NOx'])):\n", "\n", " time = d['time']\n", "\n", " units = 'Mt a$^{-1}$' if gas == 'CO2' else 'kt a$^{-1}$'\n", " factor = 1e3 if gas == 'NOx' else 1.0\n", "\n", " values = factor * d[f'{gas}_estimates'].sel(source=source)\n", " values_std = factor * d[f'{gas}_estimates_std'].sel(source=source)\n", "\n", " #true = ddeq.smartcarb.read_true_emissions('NO2' if gas =='NOx' else gas, source)\n", " true = ddeq.smartcarb.read_true_emissions(\n", " gas,\n", " source,\n", " time=pd.date_range('2015-01-01 11:00', '2015-12-31 11:00')\n", " )\n", " true = factor * ucat.convert_mass_per_time_unit(true, 'kg/s', 'Mt/y')\n", "\n", " mask = (values > 0) & (values < 100) & (values_std < 20)\n", "\n", " fit_result, func, times, cycle, chi2 = ddeq.timeseries.fit(\n", " time.values[mask], values[mask], values_std[mask])\n", "\n", " em, em_std = func.integrate(fit_result['x'], fit_result['x_std'])\n", "\n", " axes[i].errorbar(time[mask].values, values[mask], yerr=values_std[mask],\n", " ls='', marker='o', ms=4, label='Individual estimates')\n", "\n", " axes[i].plot(true.time, true, '-', label=f'True emissions\\n({np.mean(true):.1f} {units})')\n", " axes[i].plot(times, cycle, 'k--', label=f'Estimated emissions\\n({em:.1f}$\\\\pm${em_std:.1f} {units})')\n", " axes[i].grid(True)\n", "\n", " em, em_std = func.integrate(fit_result.x, fit_result.x_std)\n", "\n", " axes[i].set_ylim(0,100)\n", " axes[i].set_title(f'({\"abcd\"[i]}) {ddeq.vis.sub_numbers(gas)} emissions of {source}')\n", "\n", " axes[i].set_ylabel(f'{ddeq.vis.sub_numbers(gas)} [{units}]')\n", "\n", " axes[i].set_xlim(pd.Timestamp(2015,1,1), pd.Timestamp(2015,12,31))\n", " axes[i].legend(ncol=1, loc=(1.02, 0.06))\n", " axes[i].xaxis.set_major_formatter(mdates.DateFormatter('%b'))\n", "\n", "axes[-1].set_xlabel('2015')\n", "\n", "plt.tight_layout()\n", "plt.savefig('CSF-timeseries.png', dpi=500)" ] }, { "cell_type": "markdown", "id": "5674cea5-3b32-4ff6-afa0-43a46a1ed5c6", "metadata": {}, "source": [ "## Divergence method\n", "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).\n", "\n", "The code below can be used for reprocessing the example. Note that it requires the SMARTCARB dataset. The processing time is about 20 minutes." ] }, { "cell_type": "code", "execution_count": null, "id": "b7181260-49a1-411e-87db-9d7a36f18acd", "metadata": {}, "outputs": [], "source": [ "# on ICOS-CP computation time is about 1 minute per month and species \n", "# (i.e. processing one year of CO2 and NO2 took about 20 minutes)\n", "reprocess_div = False\n", "\n", "if reprocess_div:\n", "\n", " # change this to point to the SMARTCARB dataset for divergence method\n", " ROOT = '/project/coco2/fileshare/WP4/SMARTCARB/'\n", " SMARTCARB_DATA_PATH = os.path.join(ROOT, 'level2')\n", " SMARTCARB_WIND_PATH = os.path.join(ROOT, 'wind_fields')\n", "\n", " sources = ddeq.misc.read_point_sources()\n", " sources = sources.sel(source=['Berlin', 'Janschwalde'])\n", "\n", " l2 = ddeq.smartcarb.Level2Dataset(SMARTCARB_DATA_PATH, 'ace', co2_noise_scenario='low', no2_noise_scenario='low',\n", " make_no2_error_cloud_dependent=False)\n", "\n", " res = ddeq.div.estimate_emissions(\n", " l2, SMARTCARB_WIND_PATH, sources=sources,\n", " pattern='SMARTCARB_winds_%Y%m%d%H.nc',\n", " wind_product='SMARTCARB', trace_gases=['CO2', 'NO2'],\n", " varnames=['CO2', 'NO2'], smooth_data=[True, False],\n", " remove_background=[True, False],\n", " start_date='2015-01-01', end_date='2015-12-24'\n", " )\n", "\n", " res.to_netcdf('SMARTCARB_DIV_example.nc')" ] }, { "cell_type": "markdown", "id": "7e32d65a-3bad-47b5-8712-04772d25d677", "metadata": {}, "source": [ "Open the example file and print estimated CO2 and NOX emissions:" ] }, { "cell_type": "code", "execution_count": null, "id": "b5b0e116-63e2-4305-ace2-fd6e0ab45eab", "metadata": {}, "outputs": [], "source": [ "res = xr.open_dataset(os.path.join(ddeq.DATA_PATH, 'SMARTCARB_DIV_example.nc'))\n", "\n", "annual = {\n", " 'CO2': {},\n", " 'NO2': {},\n", "}\n", "\n", "for source in ['Janschwalde', 'Berlin']:\n", " for gas in ['CO2', 'NO2']:\n", " em = res[f'{gas}_emissions'].sel(source=source)\n", " em_std = res[f'{gas}_emissions_precision'].sel(source=source)\n", "\n", " factor = 1.0 if gas == 'CO2' else 1.32\n", "\n", " units = 'Mt' if gas == 'CO2' else 'kt'\n", " em = factor * ucat.convert_mass_per_time_unit(em, 'kg/s', f'{units}/a')\n", " em_std = factor * ucat.convert_mass_per_time_unit(em_std, 'kg/s', f'{units}/a')\n", "\n", " print(f'{source:20s}{em:5.1f} ± {em_std:5.1f} {units} {gas} per year')\n", "\n", " annual[gas][source] = f'{em:.1f} ± {em_std:.1f} {units}'" ] }, { "cell_type": "markdown", "id": "85e382ef-f5ad-4b12-8d95-42345f7fc4d9", "metadata": {}, "source": [ "The following code plots the enhancements and fluxes of CO2 and NOX over Berlin and Jänschwalde\n", "from one year of SMARTCARB data:" ] }, { "cell_type": "code", "execution_count": null, "id": "56fc400c-6bfd-4de9-89c4-7b4b6ecf687a", "metadata": {}, "outputs": [], "source": [ "DOMAIN = ddeq.misc.Domain('BER+JAE', 12.5, 51.25, 15.3, 53.1)" ] }, { "cell_type": "code", "execution_count": null, "id": "1198888b-8b15-4cc8-9bf3-cbfa219e4e24", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(2,2, figsize=(12,6.6),\n", " subplot_kw={'projection': DOMAIN.proj})\n", "\n", "# panel a\n", "fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[0,0], admin_level=1, edgecolor='k')\n", "\n", "for source in ['Berlin', 'Janschwalde']:\n", " this = res.sel(source=source)\n", " m = ax.pcolormesh(\n", " this.longrid, this.latgrid,\n", " ucat.convert_columns(this['CO2_grid'], 'm-2', 'ppmv', molar_mass='CO2'),\n", " vmin=-0.4, vmax=0.4, transform=ccrs.PlateCarree())\n", "\n", "plt.colorbar(m, ax=ax).set_label('XCO$_2$ enhancement [ppmv]')\n", "\n", "\n", "# panel b\n", "fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[0,1], admin_level=1, edgecolor='k')\n", "\n", "for source in ['Berlin', 'Janschwalde']:\n", " this = res.sel(source=source)\n", " m = ax.pcolormesh(this.longrid, this.latgrid,\n", " ucat.convert_columns(this['NO2_grid'], 'm-2', 'umol m-2', molar_mass='NO2'),\n", " vmin=0, vmax=100, transform=ccrs.PlateCarree())\n", "\n", "plt.colorbar(m, ax=ax).set_label('NO$_2$ enhancement [µmol m$^{-2}$]')\n", "\n", "\n", "# panel c\n", "fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[1,0], admin_level=1, edgecolor='k')\n", "\n", "for source in ['Berlin', 'Janschwalde']:\n", " this = res.sel(source=source)\n", " m = ax.pcolormesh(this.longrid, this.latgrid,\n", " this['CO2_div'] / 1e3,\n", " vmin=-2, vmax=4, transform=ccrs.PlateCarree())\n", "\n", "plt.colorbar(m, ax=ax).set_label('XCO$_2$ flux [kg m$^{-2}$ s$^{-1}$]')\n", "\n", "ax.text(0.975, 0.96, f'Berlin: {annual[\"CO2\"][\"Berlin\"]}\\nJänschwalde: {annual[\"CO2\"][\"Janschwalde\"]}',\n", " ha='right', va='top', fontsize='small', transform=ax.transAxes,\n", " bbox={'ec': 'black', 'fc': 'white', 'lw': 0.15}\n", " )\n", "\n", "\n", "# panel d\n", "fig, ax, _ = ddeq.vis.create_map(DOMAIN, ax=axes[1,1], admin_level=1, edgecolor='k')\n", "\n", "for source in ['Berlin', 'Janschwalde']:\n", " this = res.sel(source=source)\n", " m = ax.pcolormesh(this.longrid, this.latgrid,\n", " this['NO2_div'],\n", " vmin=-2, vmax=4, transform=ccrs.PlateCarree())\n", "\n", "plt.colorbar(m, ax=ax).set_label('NO$_x$ flux [g m$^{-2}$ s$^{-1}$]')\n", "\n", "ax.text(0.975, 0.96, f'Berlin: {annual[\"NO2\"][\"Berlin\"]}\\nJänschwalde: {annual[\"NO2\"][\"Janschwalde\"]}',\n", " ha='right', va='top', fontsize='small', transform=ax.transAxes,\n", " bbox={'ec': 'black', 'fc': 'white', 'lw': 0.15}\n", " )\n", "\n", "for i, ax in enumerate(axes.flatten()):\n", " ddeq.vis.add_gridlines(ax, dlon=0.5, dlat=0.5)\n", " ddeq.vis.add_hot_spots(\n", " ax,\n", " sources=sources.sel(source=['Berlin', 'Janschwalde', 'Schwarze Pumpe', 'Boxberg']),\n", " do_path_effect=True, ms=3, size='small'\n", " )\n", " ax.set_title(f'({\"abcd\"[i]})')\n", "\n", "plt.tight_layout()\n", "plt.savefig('DIV-example.png', dpi=500)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.21" } }, "nbformat": 4, "nbformat_minor": 5 }