emulator_training.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "88fff337",
   "metadata": {},
   "source": [
    "# Emulation of air quality in China"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e95275b0",
   "metadata": {},
   "source": [
    "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/lukeconibear/emulator_training/blob/main/emulator_training.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8287eb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "IN_COLAB = 'google.colab' in sys.modules"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b1f33e6",
   "metadata": {},
   "source": [
    "![emulator.png](https://user-images.githubusercontent.com/19871268/128219782-a93b83a1-eafa-4e06-851e-92b4e10c256b.png)\n",
    "\n",
    "Images [1](https://www.nasa.gov/content/a-portrait-of-global-winds/), [2](https://www.aidanscannell.com/post/gaussian-process-regression/), and [3](http://6degreesoffreedom.co/circle-random-sampling/).  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7e3bf61",
   "metadata": {},
   "source": [
    "## Approach\n",
    "\n",
    "### Problem\n",
    "1. Identify\n",
    "    - What is the problem you're trying to solve?\n",
    "2. Inputs\n",
    "    - What inputs go into the emulator?\n",
    "3. Outputs\n",
    "    - What will the emulator predict?\n",
    "\n",
    "### Simulations\n",
    "4. Design\n",
    "    - Determine what the simulations will be (taking account of the inputs above).  \n",
    "5. Run\n",
    "    - Perform these simulations (including a control run).  \n",
    "6. Evaluate\n",
    "    - Compare the control simulation to measurements to ensure accurately represents the outputs above.\n",
    "\n",
    "### Emulators\n",
    "7. Design\n",
    "    - Optimise the design of the emulators.  \n",
    "8. Training\n",
    "    - Train the emulators on the simulation data.\n",
    "9. Evaluate\n",
    "    - Test the emulators on the held-out test data.\n",
    "10. Predict\n",
    "    - Use the emulators for new predictions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42e0f3e8",
   "metadata": {},
   "source": [
    "## Case study"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bde97c9",
   "metadata": {},
   "source": [
    "### 1. Problem: Identify"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db9a898b",
   "metadata": {},
   "source": [
    "Predict air quality from emission changes.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbdada3c",
   "metadata": {},
   "source": [
    "### 2. Problem: Inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85faae04",
   "metadata": {},
   "source": [
    "5 inputs, key anthropogenic emission sectors:  \n",
    "- Residential (RES)\n",
    "- Industrial (IND)\n",
    "- Land transport (TRA)\n",
    "- Agriculture (AGR)\n",
    "- Power generation (ENE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f93f4aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "sectors = ['RES', 'IND', 'TRA', 'AGR', 'ENE']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d1068a4",
   "metadata": {},
   "source": [
    "### 3. Problem: Outputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a2a040d",
   "metadata": {},
   "source": [
    "1 output: PM$_{2.5}$ concentrations.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a222b69",
   "metadata": {},
   "outputs": [],
   "source": [
    "output = 'PM2_5_DRY'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5a714ec",
   "metadata": {},
   "source": [
    "### 4. Simulations: Design"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "006b3bde",
   "metadata": {},
   "source": [
    "Maximin [Latin-hypercube design](https://en.wikipedia.org/wiki/Latin_hypercube_sampling) using [pyDOE's](https://pythonhosted.org/pyDOE/index.html) [`lhs`](https://pythonhosted.org/pyDOE/randomized.html#latin-hypercube-lhs) method.\n",
    "- For training, the number of simulation runs required is 10 times the number of inputs to adequately cover the parameter space (e.g., for 5 inputs, you'll need 50 runs).  \n",
    "- For testing, have a separate set of simulations. A smaller number will suffice (e.g., for 5 inputs, 5 test simulations).  \n",
    "- The Latin-hypercube design enables a near-random sample of values across the whole parameter space.  \n",
    "- Maxi-min: Decides on the maximum and minimum values for the inputs, defaults to 1.0 and 0.0.  \n",
    "- Choose number of iterations for the sampling (100,000 is common, though using 1,000 here is the example to run faster).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdf1f7a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "if IN_COLAB:\n",
    "    !pip install --upgrade pyDOE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe47564a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pyDOE import lhs\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec9682a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_inputs = 5\n",
    "iterations = 1_000 # use 100_000 in actual work\n",
    "\n",
    "# create latin hypercube design, sampled between 0 (min) and 1 (max)\n",
    "design_train = lhs(\n",
    "    number_inputs, \n",
    "    samples=number_inputs * 10, \n",
    "    criterion='maximin', \n",
    "    iterations=iterations)\n",
    "\n",
    "design_test = lhs(\n",
    "    number_inputs, \n",
    "    samples=number_inputs, \n",
    "    criterion='maximin', \n",
    "    iterations=iterations)\n",
    "\n",
    "# change scaling of max/min if appropriate\n",
    "max_value = 1.5 # scale inputs up to 150% of baseline\n",
    "\n",
    "design_train = design_train * max_value\n",
    "design_test = design_test * max_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac6681d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert numpy arrays to pandas dataframes\n",
    "inputs_train = pd.DataFrame(data=design_train, columns=sectors)\n",
    "inputs_test = pd.DataFrame(data=design_test, columns=sectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f67cb47",
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs_train.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa41c025",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the pairs as a histogram to check random sample over parameter space\n",
    "sns.pairplot(inputs_train, diag_kind='hist')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c56ab2b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "if IN_COLAB:\n",
    "    !wget https://raw.githubusercontent.com/lukeconibear/emulator_training/main/latin_hypercube_inputs_train.csv\n",
    "    !wget https://raw.githubusercontent.com/lukeconibear/emulator_training/main/latin_hypercube_inputs_test.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0557c15e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in actual ones used\n",
    "inputs_train = pd.read_csv('latin_hypercube_inputs_train.csv')\n",
    "inputs_test = pd.read_csv('latin_hypercube_inputs_test.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d856da05",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to arrays\n",
    "X_train, X_test = inputs_train.values, inputs_test.values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "137d6e46",
   "metadata": {},
   "source": [
    "### 5. Simulations: Run"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91bc81b8",
   "metadata": {},
   "source": [
    "Use a complex air quality model ([WRFChem](https://wrfchem-leeds.github.io/WRFotron/)).  \n",
    "\n",
    "- 1 year for the control simulation.  \n",
    "\n",
    "Then step 6, before:\n",
    "\n",
    "- 50 years of training simulations.  \n",
    "- 5 years of testing simulations. \n",
    "\n",
    "*Warning: this takes a long time, even on a [high-performance computer](https://arcdocs.leeds.ac.uk/welcome.html).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2793f0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if IN_COLAB:\n",
    "    !wget https://raw.githubusercontent.com/lukeconibear/emulator_training/main/df_train.csv\n",
    "    !wget https://raw.githubusercontent.com/lukeconibear/emulator_training/main/df_test.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea83b7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bringing over data from my WRFChem runs \n",
    "df_train = pd.read_csv('df_train.csv', index_col='Unnamed: 0')\n",
    "df_test = pd.read_csv('df_test.csv', index_col='Unnamed: 0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "792641c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "042e3ff1",
   "metadata": {},
   "outputs": [],
   "source": [
    "gridcells = df_train[['lat', 'lon']].drop_duplicates().values.tolist()\n",
    "\n",
    "print(f'There are {len(gridcells)} grid cells in the domain.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec6ffc32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick a random grid cell (no reproducible seed, lucky dip!)\n",
    "import random\n",
    "\n",
    "random_gridcell_lat, random_gridcell_lon = random.choice(gridcells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08aed9e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets see what you picked ...\n",
    "print(f'Lat = {random_gridcell_lat}')\n",
    "print(f'Lon = {random_gridcell_lon}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bd46ee2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# and what was the output for this random grid cell over all training simulations\n",
    "\n",
    "random_gridcell_y_train = df_train.loc[df_train.lat == random_gridcell_lat].loc[df_train.lon == random_gridcell_lon][[output]].values\n",
    "random_gridcell_y_test = df_test.loc[df_test.lat == random_gridcell_lat].loc[df_test.lon == random_gridcell_lon][[output]].values\n",
    "\n",
    "random_gridcell_y_train"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e40d078d",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 6. Simulations: Evaluate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58423ff7",
   "metadata": {},
   "source": [
    "Evaluate it against measurements to ensure it accurately predicts outputs.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dd336e6",
   "metadata": {},
   "source": [
    "### 7. Emulators: Design"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21c5c20c",
   "metadata": {},
   "source": [
    "The emulators are based around Gaussian process models.\n",
    "___\n",
    "\n",
    "#### Explanations and interactive material\n",
    "\n",
    "- They are explained very well in the book by [Rasmussen & Williams (2006)](http://www.gaussianprocess.org/gpml/chapters/RW.pdf).\n",
    "- Geoff Hinton explains them briefly in this short [video](https://youtu.be/i0cKa0di_lo).\n",
    "- There are many nice visual explorations of how they work, such as [here](https://distill.pub/2019/visual-exploration-gaussian-processes/), [here](http://www.infinitecuriosity.org/vizgp/), and [here](https://infallible-thompson-49de36.netlify.app/#section-5).\n",
    "- Many previous studies have used Gaussian Processes, though some have used also neural networks and other machine learning methods. [Watson-Parris (2021)](https://royalsocietypublishing.org/doi/10.1098/rsta.2020.0098) provides more information on this.\n",
    "___\n",
    "\n",
    "#### Library options\n",
    "\n",
    "There is an open-source Python library, [Earth System Emulator (ESEm)](https://esem.readthedocs.io/en/latest/), which is an excellent tool for emulating geophysical datasets including Earth System Models. This has a range of emulators, including Gaussian processes and neural networks set up for you. I highly recommend getting involved (contributing) with this project, if you're interested in emulators for the Earth sciences.\n",
    "\n",
    "There are many other options including [scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html) (used here), [GPflow](https://gpflow.readthedocs.io/en/master/), and [PyMC3](https://docs.pymc.io/en/v3/).  \n",
    "___\n",
    "\n",
    "#### Optimisation\n",
    "\n",
    "We optimised the design and hyperparameters of the model using genetic programming (automatic machine learning tool, [TPOT](http://epistasislab.github.io/tpot/)).  \n",
    "\n",
    "Our final design included:\n",
    "- Gaussian process regressors are implemented in [scikit-learn](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html).  \n",
    "- Preprocessing the inputs with a power transform ([Yeo-Johnson](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PowerTransformer.html)) to make them more Gaussian-like.  \n",
    "- We used a [Matern](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.Matern.html?highlight=matern#sklearn.gaussian_process.kernels.Matern) kernel.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df263aa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import PowerTransformer\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "from sklearn.gaussian_process.kernels import Matern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f4a0157",
   "metadata": {},
   "outputs": [],
   "source": [
    "emulator = make_pipeline(\n",
    "    PowerTransformer(),\n",
    "    GaussianProcessRegressor(\n",
    "        kernel=Matern(length_scale=3.4, nu=2.5), \n",
    "        n_restarts_optimizer=240, # to help avoid local maximums\n",
    "        normalize_y=False)\n",
    ")\n",
    "emulator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a784ea3",
   "metadata": {},
   "source": [
    "### 8. Emulators: Training"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e605aaee",
   "metadata": {},
   "source": [
    "The emulators were trained on the 50 years of simulator data (in parallel on a high-performance computer).  \n",
    "Example below for the single random grid cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd9b213a",
   "metadata": {},
   "outputs": [],
   "source": [
    "emulator.fit(X_train, random_gridcell_y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "417e40b5",
   "metadata": {},
   "source": [
    "### 9. Emulators: Evaluate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "416df1b3",
   "metadata": {},
   "source": [
    "Predict unseen test data and evaluate against simulated values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac545b05",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4350cab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def root_mean_squared_error(y_true, y_pred):\n",
    "    return np.round(np.sqrt(np.abs(mean_squared_error(y_true, y_pred))), decimals=4)\n",
    "\n",
    "def coefficient_of_determination(y_true, y_pred):\n",
    "    return np.round(r2_score(y_true, y_pred), decimals=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b444c63",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_scatter_plot(lat, lon, output, y_test, y_pred, rmse_test, r2_test):\n",
    "    # setup the plot\n",
    "    fig = plt.figure(1, figsize=(5, 5))\n",
    "    gs = gridspec.GridSpec(1, 1)\n",
    "    ax = fig.add_subplot(gs[0])\n",
    "    \n",
    "    # make it look nice\n",
    "    ax.set_facecolor('whitesmoke')\n",
    "    plt.title(f'Random grid cell\\n(lat = {lat}, lon = {lon})', fontsize=14)\n",
    "    limit = int(np.ceil(np.nanmax(random_gridcell_y_pred) / 5)) * 5\n",
    "    plt.xlim([0, limit])\n",
    "    plt.ylim([0, limit])\n",
    "    ticks = np.linspace(0, 20, limit // int(limit / 5))\n",
    "    plt.xticks(ticks, fontsize=14)\n",
    "    plt.yticks(ticks, fontsize=14)\n",
    "    label = 'annual-mean PM$_{2.5}$\\nconcentrations (${\\mu}g$ $m^{-3}$)'\n",
    "    plt.xlabel('Simulator, ' + label, fontsize=14)\n",
    "    plt.ylabel('Emulator, ' + label, fontsize=14)\n",
    "    x = np.arange(2 * np.ceil(limit))\n",
    "    plt.plot(x, x, color='grey', linestyle='--')\n",
    "    plt.plot(x, 0.5 * x, color='grey', linestyle='--')\n",
    "    plt.plot(x, 2 * x, color='grey', linestyle='--')\n",
    "    text = \"R$^2$ = \" + str(r2_test) + \"\\nRMSE = \" + str(rmse_test)\n",
    "    at = matplotlib.offsetbox.AnchoredText(text, prop=dict(size=14), frameon=True, loc='upper left')\n",
    "    at.patch.set_boxstyle(\"round,pad=0.,rounding_size=0.2\")\n",
    "    ax.add_artist(at)\n",
    "    \n",
    "    # show the data\n",
    "    plt.scatter(np.vstack(y_test), np.vstack(y_pred), s=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f487e5f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "random_gridcell_y_pred = emulator.predict(X_test)\n",
    "random_gridcell_y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4d44da0",
   "metadata": {},
   "outputs": [],
   "source": [
    "random_gridcell_rmse_test = root_mean_squared_error(random_gridcell_y_test, random_gridcell_y_pred)\n",
    "random_gridcell_r2_test = coefficient_of_determination(random_gridcell_y_test, random_gridcell_y_pred)\n",
    "\n",
    "print(f'The random grid cell has a RSME value of: {random_gridcell_rmse_test} ugm-3')\n",
    "print(f'The random grid cell has a R^2 value of:  {random_gridcell_r2_test}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c59143d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# example for the random grid cell\n",
    "create_scatter_plot(\n",
    "    random_gridcell_lat,\n",
    "    random_gridcell_lon,\n",
    "    output, \n",
    "    random_gridcell_y_test, \n",
    "    random_gridcell_y_pred, \n",
    "    random_gridcell_rmse_test,\n",
    "    random_gridcell_r2_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31828d7e",
   "metadata": {},
   "source": [
    "We developed 1 emulator per grid cell from the simulator (15,278 grid cells) to capture the spatial distribution of the output."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca867027",
   "metadata": {},
   "source": [
    "### 10. Emulators: Predict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06c1e584",
   "metadata": {},
   "source": [
    "Create your own custom emission configuration and use your individual emulator to predict what the PM$_{2.5}$ concentrations would be. Remember, this is still just for a single grid cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6602774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select values between 0.0 and 1.5\n",
    "fraction_of_residential_emissions      = 0.6\n",
    "fraction_of_industrial_emissions       = 1.4\n",
    "fraction_of_land_transport_emissions   = 1.3\n",
    "fraction_of_agricultural_emissions     = 1.2\n",
    "fraction_of_power_generation_emissions = 0.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc7e539a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# placing these into a stacked column array for the emulator\n",
    "custom_emission_configuration = np.array([[\n",
    "    fraction_of_residential_emissions,\n",
    "    fraction_of_industrial_emissions,\n",
    "    fraction_of_land_transport_emissions,\n",
    "    fraction_of_agricultural_emissions,\n",
    "    fraction_of_power_generation_emissions]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9477b6b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now predict the output for these inputs\n",
    "prediction = emulator.predict(custom_emission_configuration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56f2690e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is just some setup for the table and plot below\n",
    "from tabulate import tabulate\n",
    "\n",
    "table = tabulate(\n",
    "    [\n",
    "        ['Residential', int((fraction_of_residential_emissions * 100) - 100)],\n",
    "        ['Industrial', int((fraction_of_industrial_emissions * 100) - 100)],\n",
    "        ['Land transport', int((fraction_of_land_transport_emissions * 100) - 100)],\n",
    "        ['Agricultural', int((fraction_of_agricultural_emissions * 100) - 100)],\n",
    "        ['Power generation', int((fraction_of_power_generation_emissions * 100) - 100)],\n",
    "    ],\n",
    "    headers=['Emission sector', 'Fractional reduction (%)'])\n",
    "\n",
    "# ---\n",
    "import folium\n",
    "\n",
    "background_map = folium.Map(\n",
    "    location=[random_gridcell_lat, random_gridcell_lon], \n",
    "    zoom_start=5, \n",
    "    tiles=\"Stamen Terrain\")\n",
    "\n",
    "tooltip = \"Click me!\"\n",
    "\n",
    "folium.Marker(\n",
    "    [random_gridcell_lat, random_gridcell_lon],\n",
    "    popup=f\"<b>Emulator</b> {float(prediction):.1f} \\u03BCg/m\\u00b3\",\n",
    "    tooltip=tooltip\n",
    ").add_to(background_map);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "842445a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'You choose a custom emission configuration where:')\n",
    "print()\n",
    "print(table)\n",
    "print()\n",
    "print(f'The predicted PM2.5 concentration would be {float(prediction):.1f}' + ' \\u03BCg/m\\u00b3')\n",
    "print()\n",
    "background_map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71c992c9",
   "metadata": {},
   "source": [
    "All 15,278 emulators were then used to predict air quality concentrations for all emission configurations within a 0−150% matrix of emission scaling factors at 20% increments (32,768 emission configurations).  \n",
    "\n",
    "For example, one emission configuration could be (*emulator inputs*):  \n",
    "\n",
    "| Sector           | Fractional emissions              |\n",
    "| ---              | ---                               |\n",
    "| Residential      | 0.6 (-40% from the 2015 baseline) | \n",
    "| Industrial       | 0.4 (-60% from the 2015 baseline) | \n",
    "| Land transport   | 0.8 (-20% from the 2015 baseline) | \n",
    "| Agricultural     | 0.8 (-20% from the 2015 baseline) | \n",
    "| Power generation | 0.8 (-20% from the 2015 baseline) | \n",
    "\n",
    "This result of this emission configuration should be (*emulator outputs*):\n",
    "- PM$_{2.5}$ exposure in China of 34.3 ${\\mu}g$ $m^{-3}$ (-28%).  \n",
    "    - Reaching the National Air Quality Target (NAQT) of 35 ${\\mu}g$ $m^{-3}$.  \n",
    "- 361,900 avoided premature deaths per year (-17%).\n",
    "    - 1,781,800 remaining premature deaths per year.\n",
    "\n",
    "Play around with the interactive plot to see what air quality and public health benefits you can find.  \n",
    "\n",
    "*The code for the interactive plot is [here](https://github.com/lukeconibear/emulator_training/blob/main/interactive_plot.py).* \n",
    "\n",
    "**Questions?**  \n",
    "- Removing emissions from which 1 sector would attain the NAQT?\n",
    "- Removing emissions from which 2 sectors would attain the World Health Organization (WHO) Interim Target 2 (25 ${\\mu}g$ $m^{-3}$)?\n",
    "- Does removing emissions from all of these 5 sectors attain the WHO Annual Guideline (5 ${\\mu}g$ $m^{-3}$)?\n",
    "    - Why?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be91b5cf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "\n",
    "IFrame(\"https://emulator.pythonanywhere.com/\", width=800, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b313fbb-98f8-4157-b397-a16d6cb65820",
   "metadata": {},
   "source": [
    "## Possible next steps\n",
    "- Play around with:\n",
    "    - The simple emulator in this notebook.\n",
    "    - The interactive guides in these articles: [here](https://distill.pub/2019/visual-exploration-gaussian-processes/), [here](http://www.infinitecuriosity.org/vizgp/), and [here](https://infallible-thompson-49de36.netlify.app).  \n",
    "- Try out the [Earth System Emulator (ESEm)](https://esem.readthedocs.io/en/latest/).  \n",
    "- Check out some climate AI communities e.g., [Climate Informatics](http://www.climateinformatics.org/), [ClimateChange AI](https://www.climatechange.ai/)\n",
    "- Read the Gaussian Process book by [Rasmussen & Williams (2006)](http://www.gaussianprocess.org/gpml/chapters/RW.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abafba41",
   "metadata": {},
   "source": [
    "## Further information\n",
    "\n",
    "### Short-term air quality prediction.\n",
    "- [Code](https://github.com/lukeconibear/emulator).\n",
    "- Papers:\n",
    "    - *Conibear, L. Reddington, C. L., Silver, B. J., Chen, Y., Knote, C., Arnold, S. R., Spracklen, D. V. (2021). Statistical emulation of winter ambient fine particulate matter concentrations from emission changes in China, GeoHealth, [DOI](https://doi.org/10.1029/2021GH000391).*\n",
    "\n",
    "### Long-term air quality and public health prediction\n",
    "- [Code](https://github.com/lukeconibear/emulator_annual).\n",
    "- Papers:\n",
    "    - *Conibear, L. Reddington, C. L., Silver, B. J., Chen, Y., Knote, C., Arnold, S. R., Spracklen, D. V. (2021, in prep.). Machine Learning Emulation of Long−Term Air Pollution Exposure and Disease Burden from Emission Changes in China.*\n",
    "    - *Conibear, L. Reddington, C. L., Silver, B. J., Arnold, S. R., Turnock, S. T., Klimont, Z., Spracklen, D. V. (2021, in prep.). The Contribution of Emission Sources to the Future Air Pollution Disease Burden in China.*\n",
    "    - *Conibear, L. Reddington, C. L., Silver, B. J., Chen, Y., Knote, C., Arnold, S. R., Spracklen, D. V. (2021, in prep.). Emission Sector Impacts on Air Quality and Public Health in China from 2010−2020.*\n",
    "    \n",
    "### Health impact assessment\n",
    "Used the Global Exposure Mortality Model (GEMM) from [Burnett et al., (2018)](https://doi.org/10.1073/pnas.1803222115). See, our papers for more detail on this method, and the code [here](https://github.com/lukeconibear/health_impact_assessment)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "website",
   "language": "python",
   "name": "website"
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}