{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Tutorial: Data-driven emission quantification of hot spots\n",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "\n",
    "# import and setup matplotlib\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "\n",
    "from ddeq import DATA_PATH\n",
    "from ddeq.smartcarb import DOMAIN\n",
    "import ddeq\n",
    "\n",
    "CRS = ccrs.UTM(32)\n",
    "WGS84 = ccrs.PlateCarree()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of point sources\n",
    "sources = ddeq.misc.read_point_sources()\n",
    "sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic satellite observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = os.path.join(DATA_PATH,'Sentinel_7_CO2_2015042311_o1670_l0483.nc')\n",
    "data_level2 = ddeq.smartcarb.read_level2(filename, co2_noise_scenario='medium',\n",
    "                                         no2_noise_scenario='high')\n",
    "data_level2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.visualize(data_level2, data_level2.NO2, gas='NO2', domain=DOMAIN,\n",
    "                         sources=sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to read and visualize the (vertically integrated) COSMO-GHG fields:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = pd.Timestamp(data_level2.time.to_pandas())\n",
    "filename = os.path.join(DATA_PATH, 'cosmo_2d_2015042311.nc')\n",
    "\n",
    "data_cosmo = ddeq.smartcarb.read_cosmo(filename, 'CO2')\n",
    "\n",
    "fig = ddeq.vis.make_field_map(data_cosmo, trace_gas='CO2', domain=DOMAIN,\n",
    "                              vmin=404, vmax=408, alpha=data_cosmo['CLCT'],\n",
    "                              border=50.0, label='XCO$_2$ [ppm]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise \n",
    "* 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.\n",
    "* 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).\n",
    "* 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.\n",
    "* 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A solution can be found in the `example-hakkarainen-2022-fig-01.ipynb` file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wind fields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "winds = ddeq.wind.read_at_sources(time, sources, radius=0.05,\n",
    "                                  product='SMARTCARB', data_path=DATA_PATH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to show the wind direction at the source location providing `winds` as an argument to `ddeq.vis.visualize`: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.visualize(data_level2, data_level2.NO2, gas='NO2', domain=DOMAIN,\n",
    "                         sources=sources, winds=winds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plume detection algorithm\n",
    "Plumes are regions how satellite pixels where CO2/NO2 values are significantly enhanced above the background\n",
    "\\begin{equation}\n",
    "SNR = \\frac{X - X_\\mathrm{bg}}{\\sqrt{\\sigma_\\mathrm{random}^2 + \\sigma_\\mathrm{sys}^2}} \\geq z_\\mathrm{thr}\n",
    "\\end{equation}\n",
    "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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = os.path.join(DATA_PATH, 'Sentinel_7_CO2_2015042311_o1670_l0483.nc')\n",
    "\n",
    "data = ddeq.smartcarb.read_level2(filename, co2_noise_scenario='low',\n",
    "                                  no2_noise_scenario='high',\n",
    "                                  co_noise_scenario='low',\n",
    "                                  only_observations=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = ddeq.dplume.detect_plumes(data, sources,\n",
    "                                 variable='NO2', variable_std='NO2_std',\n",
    "                                 filter_type='gaussian', filter_size=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plume detection can be visualized using the `ddeq.vis.visualize` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.visualize(data, data['NO2'], gas='NO2', domain=DOMAIN,\n",
    "                         winds=winds, do_zoom=False, show_clouds=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise \n",
    "* Detect emission plumes Berlin and Jänschwalde using low-noise CO2 observations.\n",
    "* Increase the size of the Gaussian filter to increase number of detected pixels.\n",
    "* Visualize the results using `ddeq.vis.visualize`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Center lines and polygons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = os.path.join(DATA_PATH, 'Sentinel_7_CO2_2015042311_o1670_l0483.nc')\n",
    "\n",
    "data = ddeq.smartcarb.read_level2(\n",
    "    filename, co2_noise_scenario='low', no2_noise_scenario='high',\n",
    "    co_noise_scenario='low', only_observations=False\n",
    ")\n",
    "\n",
    "data = ddeq.dplume.detect_plumes(\n",
    "    data, sources.sel(source=['Berlin', 'Janschwalde']),\n",
    "    variable='NO2', variable_std='NO2_std',\n",
    "    filter_type='gaussian', filter_size=0.5\n",
    ")\n",
    "\n",
    "# setting `do_zoom` to True will zoom on detected plumes and `draw_gridlines` \n",
    "# will add lines\n",
    "ddeq.vis.visualize(\n",
    "    data, 'NO2', gas='NO2', domain=DOMAIN, winds=winds, do_zoom=True,\n",
    "    show_clouds=True, draw_gridlines=True\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data, curves = ddeq.plume_coords.compute_plume_line_and_coords(\n",
    "    data, crs=CRS, radius=25e3, plume_area='area'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code shows the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddeq.vis.show_detected_plumes(\n",
    "    data, curves, data['NO2'], gas='NO2', domain=DOMAIN, winds=winds,\n",
    "    do_zoom=True, crs=CRS\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare emission quantification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = ddeq.emissions.prepare_data(data, 'CO2')\n",
    "data = ddeq.emissions.prepare_data(data, 'NO2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to visualize different variables using the `variable` parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddeq.vis.visualize(\n",
    "    data, 'CO2_minus_estimated_background_mass', gas='CO2', domain=DOMAIN,\n",
    "    winds=None, do_zoom=True, show_clouds=True, draw_gridlines=True,\n",
    "    vmin=-20e-3, vmax=40e-3,\n",
    "    label='CO$_2$ enhancement [kg m$^{-2}$]'\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-sectional flux method\n",
    "The following code estimated CO$_2$ and NO$_x$ emissions for a point source (Jänschwalde) and a city (Berlin):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-sectional flux method for point source"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "winds = ddeq.wind.read_at_sources(data.time, data, product='ERA5', era5_prefix='SMARTCARB',\n",
    "                                  data_path=DATA_PATH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cross-sectional flux (csf) method is performed by the following function.\n",
    "\n",
    "Note that `f_model` gives the factor for converting NO2 to NOx line densities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = ddeq.csf.estimate_emissions(data, winds, sources, curves,\n",
    "                                      method='gauss', gases=['CO2', 'NO2'], crs=CRS,\n",
    "                                      f_model=1.32\n",
    "                                     )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results can be visualized with the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with xr.set_options(keep_attrs=True):\n",
    "    fig = ddeq.vis.plot_csf_result(\n",
    "        ['CO2', 'NO2'],\n",
    "        data, winds, results,\n",
    "        source='Janschwalde',\n",
    "        curves=curves,\n",
    "        domain=DOMAIN, crs=CRS\n",
    "    )\n",
    "\n",
    "plt.savefig('ddeq-example.png',  bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-sectional flux method for a city\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.plot_csf_result(\n",
    "    ['CO2', 'NO2'],\n",
    "    data, winds, results,\n",
    "    source='Berlin',\n",
    "    curves=curves, \n",
    "    domain=DOMAIN, crs=CRS,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Write code to quantify the CO2 and NOx emissions of other point sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrated mass enhancement\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "winds = ddeq.wind.read_at_sources(data.time, data, product='SMARTCARB',\n",
    "                                  data_path=DATA_PATH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = ddeq.ime.estimate_emissions(data, winds, sources, curves, gas='CO2')\n",
    "\n",
    "print(' ' * 10, '\\tEstimate\\tTrue')\n",
    "\n",
    "for name, source in sources.groupby('source'):\n",
    "    if name in data.source:\n",
    "        Q = results['CO2_estimated_emissions'].sel(source=name).values\n",
    "        Q_true = ddeq.smartcarb.read_true_emissions(\n",
    "            pd.Timestamp(data.time.values), 'CO2', name\n",
    "        ).mean()\n",
    "        print(f'{name:10s}\\t{ddeq.misc.kgs_to_Mtyr(Q):.1f} Mt/a\\t{Q_true:.1f} Mt/a')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = ddeq.ime.estimate_emissions(data, winds, sources, curves,\n",
    "                                      gas='NO2', decay_time=4*60**2)\n",
    "\n",
    "results = ddeq.emissions.convert_NO2_to_NOx_emissions(results, f=1.32)\n",
    "\n",
    "print(' ' * 10, '\\tEstimate\\tTrue')\n",
    "\n",
    "for name, source in sources.groupby('source'):\n",
    "    if name in data.source:\n",
    "        Q = results['NOx_estimated_emissions'].sel(source=name).values\n",
    "        Q_true = ddeq.smartcarb.read_true_emissions(\n",
    "            pd.Timestamp(data.time.values), 'NO2', name\n",
    "        ).mean()\n",
    "        print(f'{name:10s}\\t{ddeq.misc.kgs_to_Mtyr(Q*1000):.1f} kt/a\\t{Q_true:.1f} kt/a')"
   ]
  }
 ],
 "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.10.6"
  },
  "toc-autonumbering": true,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}