{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d8002a9b-edff-4381-8dfd-1568156ae3fa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Example: TROPOMI NO$_2$ plume\n",
    "An example demonstrating the application of the CSF and Gaussian plume inversion for estimating NOx emissions from the Matimba/Medupi power plant using TROPOMI NO2 images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0a070ea-f8aa-431e-8d0d-2c73dcb6137a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import cartopy.crs as ccrs\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",
    "sources = ddeq.misc.read_point_sources()\n",
    "filename = os.path.join(ddeq.DATA_PATH, 'Matimba_S5P_RPRO_L2__NO2____20210725T110715.nc')\n",
    "\n",
    "DOMAIN = ddeq.misc.Domain('Matimba', 25, -25.2, 29, -22.9)\n",
    "CRS = ccrs.epsg(22293)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2422f026-bbc8-45d3-a653-0cfdcb6c6799",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = xr.open_dataset(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "918df6dc-c7d3-4362-91c7-faaf75859c04",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.visualize(\n",
    "    data, 1e6 * data['NO2'],\n",
    "    sources=sources,\n",
    "    vmin=0,\n",
    "    vmax=300,\n",
    "    label=\"NO$_2$ vertical column density [µmol m$^{-2}$]\",\n",
    "    units='umol m-2'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e47c234b-2bd6-48c0-b1f3-b6ebe5ddfe29",
   "metadata": {},
   "source": [
    "## Cross-sectional flux method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92cfce6c-e759-40ae-8b9a-2dec972bab85",
   "metadata": {},
   "outputs": [],
   "source": [
    "var_sys = (ucat.convert_columns(0.5e15, 'cm-2', 'mol m-2',\n",
    "                                molar_mass='NO2'))**2\n",
    "\n",
    "data = ddeq.dplume.detect_plumes(data, sources,\n",
    "                                 variable='NO2',\n",
    "                                 variable_std='NO2_std',\n",
    "                                 var_sys=var_sys,\n",
    "                                 filter_type='gaussian',\n",
    "                                 filter_size=0.5)\n",
    "\n",
    "data, curves = ddeq.plume_coords.compute_plume_line_and_coords(\n",
    "    data, crs=CRS, radius=25e3, plume_area='area'\n",
    ")\n",
    "\n",
    "data = ddeq.emissions.prepare_data(data, 'NO2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8b714d6-7df4-40e8-a057-a5d739ba48c8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "winds = ddeq.wind.read_at_sources(data.time,\n",
    "                                  sources.sel(source=['Matimba']),\n",
    "                                  timesteps=12,\n",
    "                                  era5_prefix='Matimba',\n",
    "                                  data_path=ddeq.DATA_PATH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f86c5ce4-9c91-491f-a5cd-b64db01b3a9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_csf = ddeq.csf.estimate_emissions(data, winds, sources, curves, 'NO2',\n",
    "                                          f_model=2.24, crs=CRS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59154b37-ac87-49cd-a2ed-f7d33b151d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.show_detected_plumes(\n",
    "    data, curves, 1e6 * data['NO2'], gas='NO2', ld=results_csf,\n",
    "    vmin=0,\n",
    "    vmax=300,\n",
    "    label=\"NO$_2$ columns [µmol m$^{-2}$]\",\n",
    "    winds=winds,\n",
    "    domain=DOMAIN,\n",
    "    do_zoom=False, show_clouds=False, crs=CRS,\n",
    "    figwidth=4,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63188d98-434b-4c3b-bb68-70f0e436cb8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax  = plt.subplots(1, figsize=(8,2))\n",
    "ddeq.vis.plot_along_plume(ax, 'NO2', results_csf['Matimba'])\n",
    "ax.set_xlim(right=207)\n",
    "ax.set_ylim(0,150)\n",
    "plt.tight_layout()\n",
    "\n",
    "ax.legend().remove()\n",
    "ax.legend(loc=1, ncol=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4250c3c8-ebb8-4062-b28e-6d1799fc514c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes  = plt.subplots(1,3, figsize=(8,2), sharey=True)\n",
    "\n",
    "\n",
    "for i, ax in zip([1,9,19], axes):\n",
    "    r = results_csf['Matimba'].isel(along=i)\n",
    "    ddeq.vis.plot_across_section(r, gases=['NO2'], method='gauss', ax=ax,\n",
    "                                 legend='simple')\n",
    "\n",
    "    ax.text(-36, 19.5, '%d-%d km' % (r.xa/1e3, r.xb/1e3), ha='left', va='top')\n",
    "\n",
    "for ax in axes[1:]:\n",
    "    ax.set_ylim(-1, 20)\n",
    "    ax.set_ylabel('')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "449cd7bf-d4c4-4425-acc8-e324e72addfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddeq.vis.plot_csf_result(['NO2'], data, winds, results_csf, curves, 'Matimba',\n",
    "                         vmins=[0], vmaxs=[200e-6], crs=CRS);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "549bdeff-9320-4e62-8417-81860a753751",
   "metadata": {},
   "source": [
    "## Gaussian plume inversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22bff546-5254-4c3e-a9ff-fd60a014f08e",
   "metadata": {},
   "outputs": [],
   "source": [
    "priors = {'Matimba': {\n",
    "    'NO2': {'Q': 3.0,       # kg/s\n",
    "            'tau': 4*60**2  # seconds\n",
    "           }\n",
    "}}\n",
    "data, results_gauss = ddeq.gauss.estimate_emissions(data, winds, sources, curves,\n",
    "                                                    ['NO2'], priors=priors,\n",
    "                                                    fit_decay_times=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "159549c0-ae6d-4c4a-8318-3815809e4d69",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_gauss = ddeq.emissions.convert_NO2_to_NOx_emissions(results_gauss, f=2.38)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d63089c9-04e5-4e76-8d2f-968835e4e576",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.plot_gauss_result(data, results_gauss, ['Matimba'],\n",
    "                                 'NO2', curves, crs=CRS, vmin=0, vmax=10e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8725a432-e230-4f1b-9da5-4855325f3893",
   "metadata": {},
   "outputs": [],
   "source": [
    "data['NO2_plume_model'] = xr.DataArray(\n",
    "    ucat.convert_columns(data['NO2_plume_model_mass'], 'kg m-2', 'mol m-2',\n",
    "                         molar_mass='NO2'),\n",
    "    dims=data.NO2_plume_model_mass.dims\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf88bd10-59dd-42fe-aed2-9adfb463548c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ddeq.vis.show_detected_plumes(\n",
    "    data, curves,\n",
    "    1e6 * data['NO2_plume_model'].sel(source='Matimba'),\n",
    "    gas='NO2',\n",
    "    vmin=0,\n",
    "    vmax=300,\n",
    "    label=\"NO$_2$ column [µmol m$^{-2}$]\",\n",
    "    winds=winds,\n",
    "    domain=DOMAIN,\n",
    "    do_zoom=False, show_clouds=False, crs=CRS,\n",
    "    figwidth=4, add_polygons=False\n",
    ");\n",
    "\n",
    "xp, yp = results_gauss[f'NO2_curve'].sel(source='Matimba').T\n",
    "x, y = ddeq.misc.transform_coords(xp, yp, CRS, ccrs.PlateCarree(), use_xarray=False)\n",
    "fig.axes[0].plot(x, y, 'r', transform=ccrs.PlateCarree())"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}