la_aq_health_v1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "512b82e6",
   "metadata": {},
   "source": [
    "# Tutorial: Combining adjoint sensitivities with emissions data to calculate air quality health impacts in an interactive app\n",
    "### Example: LA County\n",
    "\n",
    "In this tutorial we will demonstrate how to combine results from the GEOS-Chem adjoint, known as sensitivities, with emisisons data from the US EPA's National Emission Inventory (NEI) to determine how different emission scenarios would impact air quality health impacts in Los Angles, CA. \n",
    "\n",
    "First, we will confirm that you have the correct files needed to perform this tutorial. You should have three sensitivity files of the form:\n",
    "\n",
    ">la_sens_no.nc <br>\n",
    "la_sens_o3.nc <br>\n",
    "la_sens_pm.nc <br>\n",
    ">\n",
    "\n",
    "Three files that scale exposure contributions to health impacts:\n",
    "\n",
    "> la_health_no.nc <br>\n",
    "la_health_o3.nc <br>\n",
    "la_health_pm.nc <br>\n",
    ">\n",
    "\n",
    "Three mask files to characterize areas of action:\n",
    "\n",
    "> STATE_05x06.csv <br>\n",
    "COUNTY_05x06.csv <br>\n",
    "CITY_05x06.csv <br>\n",
    ">\n",
    "\n",
    "And eleven emissions files each of the form:\n",
    "\n",
    "> nei_emis_onroad.nc\n",
    "\n",
    "For each of the emissions files the component of the filename following \"nei_emis_\" will indicate the sector. For all files but the masks there will be date and time information in the filenames to indicate when the file was created. These files should be stored in the same directory as this tutorial.\n",
    "\n",
    "We will present a quick example of how these files can be used in the form of an interactive applet. We will check to make sure these files are in the current directory but first we need to load in the python libraries we'll be using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5aa2f0c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for this tutorial we will use netCDF4, xarray, numpy, os, ipywidgets, and matplotlib.pyplot\n",
    "import netCDF4 as nc\n",
    "import ipywidgets as widgets\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import datetime\n",
    "from ipywidgets import interactive, interact, Layout, VBox, HBox\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "# set standard font size to 20\n",
    "plt.rcParams.update({'font.size': 20}) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b291453e",
   "metadata": {},
   "source": [
    "Let's take a look at the files in the current directory. We'll also store specific filenames used in this tutorial in variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6afc174c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a list of the files in the current directory\n",
    "path = './outputs/'\n",
    "files = os.listdir(path)\n",
    "\n",
    "# loop through these files\n",
    "for name in files: \n",
    "    if name.endswith('.nc'): \n",
    "        # save the filenames into variables for sensitivities\n",
    "        if \"_sens_PM\"   in name: pm_sens_filename = path + name\n",
    "        if \"_sens_O3\"   in name: o3_sens_filename = path + name\n",
    "        if \"_sens_NO\"   in name: no_sens_filename = path + name     \n",
    "        # for health scaling\n",
    "        if \"_health_PM\" in name: pm_health_filename = path + name\n",
    "        if \"_health_O3\" in name: o3_health_filename = path + name\n",
    "        if \"_health_NO\" in name: no_health_filename = path + name  \n",
    "        # for uncertainty\n",
    "        if \"_health_lb_PM\" in name: pm_health_lb_filename = path + name\n",
    "        if \"_health_lb_O3\" in name: o3_health_lb_filename = path + name\n",
    "        if \"_health_lb_NO\" in name: no_health_lb_filename = path + name          \n",
    "        if \"_health_ub_PM\" in name: pm_health_ub_filename = path + name\n",
    "        if \"_health_ub_O3\" in name: o3_health_ub_filename = path + name\n",
    "        if \"_health_ub_NO\" in name: no_health_ub_filename = path + name    \n",
    "            \n",
    "            \n",
    "path = './static/nei/'\n",
    "files = os.listdir(path)           \n",
    "for name in files: \n",
    "    if name.endswith('.nc'):             \n",
    "        # for emissions    \n",
    "        if \"_onroad_\" in name:\n",
    "            if \"_onroad_catx\" in name:\n",
    "                onrc_emis_filename = path + name\n",
    "            else:\n",
    "                onr_emis_filename = path + name\n",
    "        if \"_nonroad_\"    in name: nonr_emis_filename = path + name               \n",
    "        if \"_ag\"          in name: agr_emis_filename  = path + name            \n",
    "        if \"_ptnonipm\"    in name: ind_emis_filename  = path + name      \n",
    "        if \"_egu_\"        in name: egu_emis_filename  = path + name  \n",
    "            \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502f30c7",
   "metadata": {},
   "source": [
    "For the purpose of this tutorial we only consider a few emissions groups associated with transportation, agriculture, energy, and industry. Next we load three masks files that isolate emissions from specific areas allowing us to quantify the impact of emission changes at different spatial scales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "92ded715",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read masks\n",
    "MU = np.loadtxt(\"outputs/CITY_05x06.csv\", comments=\"#\", delimiter=\",\", unpack=False)\n",
    "SU = np.loadtxt(\"outputs/STATE_05x06.csv\", comments=\"#\", delimiter=\",\", unpack=False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d86adf0",
   "metadata": {},
   "source": [
    "We've loaded three masks files to quantify municipal (MU), county (CU), and state (SU) impacts. In the app we aim to consider areas of growing emissions actions so we combine the masks for larger and larger areas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "06bd7583",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The state mask is equivalent to the city + county + state. Areas that are double counted are set to one.\n",
    "SU = MU  + SU\n",
    "SU[SU > 1] = 1\n",
    "# Create a mask of ones for the whole domain\n",
    "DU = SU * 0 +1\n",
    "\n",
    "names = ['Municipal', 'State','North America']\n",
    "masks = [MU, SU, DU]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7497e73c",
   "metadata": {},
   "source": [
    "Next we read in all of the emissions, sensitivities, and health scaling data. We combine all transportation emissions (from nonroad, onroad, and onroad_catx) into a single emissions group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6a164bcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in emissions data from the onroad, onroad_catx, and nonroad files and combine them\n",
    "onr = xr.open_dataset(onr_emis_filename) + xr.open_dataset(onrc_emis_filename) + xr.open_dataset(nonr_emis_filename)\n",
    "# read in emissions data from the agriculture file\n",
    "agr = xr.open_dataset(agr_emis_filename)\n",
    "# read in emissions data from the industry file\n",
    "ind = xr.open_dataset(ind_emis_filename)\n",
    "# read in emissions data from the egu file\n",
    "egu = xr.open_dataset(egu_emis_filename)\n",
    "# read in adjoint sensitivity for PM2.5, O3, and NO2 to precursor species\n",
    "pm_sens = xr.open_dataset(pm_sens_filename)\n",
    "o3_sens = xr.open_dataset(o3_sens_filename)\n",
    "no_sens = xr.open_dataset(no_sens_filename)\n",
    "# read in health scaling files for PM2.5, O3, and NO2\n",
    "pm_health = xr.open_dataset(pm_health_filename)\n",
    "o3_health = xr.open_dataset(o3_health_filename)\n",
    "no_health = xr.open_dataset(no_health_filename)\n",
    "# read in uncertainty files for PM2.5, O3, and NO2\n",
    "pm_health_lb = xr.open_dataset(pm_health_lb_filename)\n",
    "o3_health_lb = xr.open_dataset(o3_health_lb_filename)\n",
    "no_health_lb = xr.open_dataset(no_health_lb_filename)\n",
    "pm_health_ub = xr.open_dataset(pm_health_ub_filename)\n",
    "o3_health_ub = xr.open_dataset(o3_health_ub_filename)\n",
    "no_health_ub = xr.open_dataset(no_health_ub_filename)\n",
    "# create an array of sectors emissions\n",
    "sectors = [onr,egu,agr,ind]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4689f0",
   "metadata": {},
   "source": [
    "In the next cell we calculate the health impact contributions from these four sectors for each of the three pollutants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "89f482e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize PM2.5 contributions by taking the BC sens array,\n",
    "# expand to a fourth sectoral dimension and multiply by 0\n",
    "# using xarray expand_dims not numpy\n",
    "pm_cont = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "pm_cont_lb = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "pm_cont_ub = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "# get precursor species variables loop through species\n",
    "variables = [species for species in pm_sens]\n",
    "for var in variables:\n",
    "    # erase \"_PM\" from variable name for emissions\n",
    "    variable = var.replace('_PM','') \n",
    "    # loop through sectors to calculate sectoral contribution\n",
    "    for i,sect in enumerate(sectors):\n",
    "        pm_cont[i,:,:,:] = pm_cont[i,:,:,:] + np.multiply(np.multiply(sect[variable],pm_sens[var]),pm_health[var])\n",
    "        pm_cont_lb[i,:,:,:] = pm_cont_lb[i,:,:,:] + np.multiply(np.multiply(sect[variable],pm_sens[var]),pm_health_lb[var])\n",
    "        pm_cont_ub[i,:,:,:] = pm_cont_ub[i,:,:,:] + np.multiply(np.multiply(sect[variable],pm_sens[var]),pm_health_ub[var])\n",
    "        \n",
    "# perform same calculation for O3\n",
    "o3_cont = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "o3_cont_lb = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "o3_cont_ub = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "variables = [species for species in o3_sens]\n",
    "for var in variables:\n",
    "    variable = var.replace('_O3','') \n",
    "    for i,sect in enumerate(sectors):\n",
    "        o3_cont[i,:,:,:] = o3_cont[i,:,:,:] + np.multiply(np.multiply(sect[variable],o3_sens[var]),o3_health[var])         \n",
    "        o3_cont_lb[i,:,:,:] = o3_cont_lb[i,:,:,:] + np.multiply(np.multiply(sect[variable],o3_sens[var]),o3_health_lb[var])\n",
    "        o3_cont_ub[i,:,:,:] = o3_cont_ub[i,:,:,:] + np.multiply(np.multiply(sect[variable],o3_sens[var]),o3_health_ub[var])\n",
    "\n",
    "        \n",
    "              \n",
    "# perform same calculation for NO2            \n",
    "no_cont = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "no_cont_lb = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "no_cont_ub = pm_sens.BC.expand_dims({'sectors':4}) * 0\n",
    "variables = [species for species in no_sens]\n",
    "for var in variables:\n",
    "    variable = var.replace('_NO2','') \n",
    "    for i,sect in enumerate(sectors):\n",
    "        no_cont[i,:,:,:] = no_cont[i,:,:,:] + np.multiply(np.multiply(sect[variable],no_sens[var]),no_health[var])         \n",
    "        no_cont_lb[i,:,:,:] = no_cont_lb[i,:,:,:] + np.multiply(np.multiply(sect[variable],no_sens[var]),no_health_lb[var])\n",
    "        no_cont_ub[i,:,:,:] = no_cont_ub[i,:,:,:] + np.multiply(np.multiply(sect[variable],no_sens[var]),no_health_ub[var])\n",
    "\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58ed9817",
   "metadata": {},
   "source": [
    "Next we create the four slider widgets along with the pollutant drop-down selector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e76e5ea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create slider widgets  \n",
    "a=widgets.IntSlider(min=-100, max=100, step=5, value=0,description='Transport:',layout=widgets.Layout(width='100%')) \n",
    "clear_output()\n",
    "b=widgets.IntSlider(min=-100, max=100, step=5, value=0,description='Energy:',layout=widgets.Layout(width='100%'))\n",
    "c=widgets.IntSlider(min=-100, max=100, step=5, value=0,description='Agriculture:',layout=widgets.Layout(width='100%'))\n",
    "d=widgets.IntSlider(min=-100, max=100, step=5, value=0,description='Industry:',layout=widgets.Layout(width='100%'))\n",
    "# create species selector \n",
    "species = widgets.Select(options=['PM25','O3','NO2'],description='Species:',disabled=False)\n",
    "# specify handle colors of sliders\n",
    "a.style.handle_color = 'lightgreen' \n",
    "b.style.handle_color = 'lightblue' \n",
    "c.style.handle_color = 'red' \n",
    "d.style.handle_color = 'black'  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa72ef0e",
   "metadata": {},
   "source": [
    "The function \"f\" takes inputs of the slider values from a,b,c, and d along with the species selected by the dropdown menu and calculates the health impacts associated with the scenario from the user input sliders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c2e01af5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(a, b, c, d, species): \n",
    "    # create a list of scale values from widgets\n",
    "    scales = [a,b,c,d]\n",
    "    # species specific information\n",
    "    if species == 'PM25':\n",
    "        J = pm_cont\n",
    "        J_lb = pm_cont_lb\n",
    "        J_ub = pm_cont_ub        \n",
    "        ymax = 3000\n",
    "        lab = ' Deaths'\n",
    "        longlab = 'Premature Deaths'\n",
    "        levels = np.linspace(-30,30,21)\n",
    "    if species == 'O3':\n",
    "        J = o3_cont\n",
    "        J_lb = o3_cont_lb\n",
    "        J_ub = o3_cont_ub          \n",
    "        ymax = 1000\n",
    "        lab = ' Deaths'        \n",
    "        longlab = 'Premature Deaths'        \n",
    "        levels = np.linspace(-10,10,21)\n",
    "    if species == 'NO2':\n",
    "        J = no_cont\n",
    "        J_lb = no_cont_lb\n",
    "        J_ub = no_cont_ub          \n",
    "        ymax = 35000\n",
    "        lab = ' Cases'        \n",
    "        longlab = 'New Pediatric Asthma Cases'\n",
    "        levels = np.linspace(-150,150,21)\n",
    "    # intialize and calculate contribution for pollutant health impacts\n",
    "    dJ = J[0,:,:,:].sum(axis=0)*0                \n",
    "    dJ_lb = J_lb[0,:,:,:].sum(axis=0)*0 \n",
    "    dJ_ub = J_ub[0,:,:,:].sum(axis=0)*0     \n",
    "    for i, scale in enumerate(scales):\n",
    "        dJ = dJ + J[i,:,:,:].sum(axis=0)*scale\n",
    "        dJ_lb = dJ_lb + J_lb[i,:,:,:].sum(axis=0)*scale\n",
    "        dJ_ub = dJ_ub + J_ub[i,:,:,:].sum(axis=0)*scale        \n",
    "    # intialize and calculate contributions from specific regions\n",
    "    values = [0,0,0]\n",
    "    values_lb = [0,0,0]\n",
    "    values_ub = [0,0,0]    \n",
    "    for i,mask in enumerate(masks):\n",
    "        values[i] = np.multiply(dJ,mask).sum()/100\n",
    "        values_lb[i] = np.multiply(dJ_lb,mask).sum()/100\n",
    "        values_ub[i] = np.multiply(dJ_ub,mask).sum()/100        \n",
    "    # create bar plot\n",
    "#    plt.close()\n",
    "#    fig, ax = plt.subplots(figsize=(10,20))  \n",
    "    fig = plt.figure(figsize=(10,20))\n",
    "    ax = fig.add_subplot(2, 1, 1)\n",
    "    bar1 = plt.bar(names,values,color=['cyan', 'red', 'green', 'blue', 'cyan'])    \n",
    "    error = [[0 for x in range(3)] for y in range(2)]\n",
    "    error[0][0]=values[0]-values_lb[0]\n",
    "    error[1][0]=values_ub[0]-values[0]\n",
    "    error[0][1]=values[1]-values_lb[1]\n",
    "    error[1][1]=values_ub[1]-values[1]       \n",
    "    error[0][2]=values[2]-values_lb[2]\n",
    "    error[1][2]=values_ub[2]-values[2]       \n",
    "    plt.errorbar(names,values,error,ls='none',capsize=10,ecolor=[0,0,0])\n",
    "    plt.title('Pollution Health Impacts in City')\n",
    "    plt.ylabel('Change in Annual '+longlab)\n",
    "    plt.ylim(-ymax, ymax)    \n",
    "    ax.axhline(0, color='grey', linewidth=0.8)    \n",
    "    # Add counts above the two bar graphs\n",
    "    for c,rect in enumerate(bar1):\n",
    "        height = rect.get_height()\n",
    "        if height < 0:\n",
    "            plt.text(rect.get_x() + rect.get_width() / 2.0, height, f'{height:.0f}'+lab, ha='center', va='top')            \n",
    "        else:\n",
    "            plt.text(rect.get_x() + rect.get_width() / 2.0, height, f'{height:.0f}'+lab, ha='center', va='bottom')        \n",
    "    \n",
    "\n",
    "    ax = fig.add_subplot(2, 1, 2)\n",
    "    plt.contourf(dJ/100,levels)\n",
    "    plt.set_cmap('coolwarm')\n",
    "    plt.colorbar()\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3be721d3",
   "metadata": {},
   "source": [
    "Lastly we set-up an interactive plot that takes the widget values and species and outputs the changing bar plot.\n",
    "\n",
    "In this plot, each of the bars refer to the area that is enacting the proposed emission change (represented by the slider values) with the number of health impacts occcurring in LA presented at the top of each of these bars.\n",
    "\n",
    "For example, the blue \"North America\" bar is the number of health impacts that would occur in LA if all of North America adopted the emission scenario proposed by the slider values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cf225dea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "802075a567154fecbe51397e96d61037",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(interactive(children=(IntSlider(value=0, description='Transport:', layout=Layout(width='100%'),…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create interactive plot    \n",
    "#plt.close('all')\n",
    "interactive_plot = interactive(f,a=a,b=b,c=c,d=d,species=species)\n",
    "clear_output()\n",
    "# specify size of plot by taking the output\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '1200px'\n",
    "output.layout.width = '800px'\n",
    "# display output in \"HBox\"\n",
    "HBox([interactive_plot])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "989c6aa9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9304e1a3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3f8f463",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}