makefiles_eachcity_from_mat.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "22f627cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "02920adbcb074eaf88518329139dbdf9",
       "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": [
    "#!/usr/bin/env python3\n",
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on Tue Nov  1 08:46:06 2022\n",
    "\n",
    "@author: schakrab\n",
    "\"\"\"\n",
    "\n",
    "# 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",
    "\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "# set standard font size to 20\n",
    "plt.rcParams.update({'font.size': 20}) \n",
    "  \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# get a list of the files in the current directory\n",
    "INIT     = 'HOU'#######THIS LINE is where the function for the Cities would be#####\n",
    "path = '/Users/schakrab/Desktop/Paper/input_folder/outputs/HOU/pyout/'####################CHANGE###########\n",
    "#path = '/Users/schakrab/Desktop/Paper/input_folder/outputs/matout/'\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 = '/Users/schakrab/Desktop/Paper/input_folder/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",
    "\n",
    "\n",
    "\n",
    "\n",
    "# read masks####################CHANGE###########\n",
    "MU = np.loadtxt(\"/Users/schakrab/Desktop/Paper/input_folder/outputs/HOU/pyout/CITY_05x06.csv\", comments=\"#\", delimiter=\",\", unpack=False)\n",
    "SU = np.loadtxt(\"/Users/schakrab/Desktop/Paper/input_folder/outputs/HOU/pyout/STATE_05x06.csv\", comments=\"#\", delimiter=\",\", unpack=False)\n",
    "####################CHANGE###########\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# 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]\n",
    "\n",
    "\n",
    "\n",
    "# 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",
    "#print(pm_sens)\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",
    "\n",
    "\n",
    "\n",
    "# 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",
    "print(pm_sens_filename) \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",
    "#print(np.shape(pm_cont),np.shape(sect[variable]),np.shape(pm_sens[var]),np.shape(pm_health[var]))\n",
    "# get precursor species variables loop through species\n",
    "variables = [species for species in pm_sens]\n",
    "print(variables)\n",
    "for var in variables:\n",
    "    # erase \"_PM\" from variable name for emissions\n",
    "    variable = var.replace('_PM','') \n",
    "    #print(variable)\n",
    "    # loop through sectors to calculate sectoral contribution\n",
    "    for i,sect in enumerate(sectors):\n",
    "     #   xxx=sect[variable]*pm_sens[var]\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",
    "#save into nc files \n",
    "\n",
    "\n",
    "base_dir = '/Users/schakrab/Desktop/Paper/input_folder/'\n",
    "#City Initials\n",
    "INIT     = 'HOU'####################CHANGE###########\n",
    "#City Name\n",
    "NAME     = 'Houston'####################CHANGE###########\n",
    "#State city is in\n",
    "STATE    = 'Texas'####################CHANGE###########\n",
    "SEX = 'both'\n",
    "IIPAR     = 91\n",
    "JJPAR     = 89\n",
    "dmax      = 396\n",
    "import pandas as pd\n",
    "from netCDF4 import Dataset,num2date,date2num\n",
    "import datetime\n",
    "yy = '2010';mm='12';\n",
    "date_series_min = pd.date_range(yy+'-'+mm, periods=396, freq='D') # This is DatetimeIndex, perhaps ordinary datetime objects are preferred...\n",
    "datesout = date_series_min.to_pydatetime() \n",
    "unout = 'days since 2010-12-01 00:00:00'\n",
    "#Set-Up\n",
    "########################################################################################################################################################\n",
    "#command lines### list of names of the cities and state ### Hardparse\n",
    "########################################################################################################################################################\n",
    "#Directories\n",
    "\n",
    "out_dir  = base_dir+'outputs/'+INIT+'/'\n",
    "fn =out_dir+INIT+'_PM_NO_O3.nc'\n",
    "ds = nc.Dataset(fn, 'w', format='NETCDF4')\n",
    "    \n",
    "    # % define dimensions\n",
    "time = ds.createDimension('time', dmax)\n",
    "lat = ds.createDimension('lat', JJPAR)\n",
    "lon = ds.createDimension('lon', IIPAR)\n",
    "sector =ds.createDimension('sector', 4)  \n",
    "\n",
    "\n",
    "# timevar = ds.createVariable('time','float32',('time'));timevar.setncattr('units',unout);\n",
    "# lats = ds.createVariable('lat', 'f4', ('lat',))\n",
    "# lons = ds.createVariable('lon', 'f4', ('lon',))\n",
    "# Sectorss = ds.createVariable('sector', 'f4', ('Sector',))\n",
    "\n",
    "#timevar[:]= date2num(datesout,unout);\n",
    "value = ds.createVariable('pm_cont', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=pm_cont\n",
    "value = ds.createVariable('pm_cont_lb', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=pm_cont_lb\n",
    "value = ds.createVariable('pm_cont_ub', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=pm_cont_ub\n",
    "value = ds.createVariable('o3_cont', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=o3_cont\n",
    "value = ds.createVariable('o3_cont_lb', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=o3_cont_lb\n",
    "value = ds.createVariable('o3_cont_ub', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=o3_cont_ub\n",
    "value = ds.createVariable('no_cont', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=no_cont\n",
    "value = ds.createVariable('no_cont_lb', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=no_cont_lb\n",
    "value = ds.createVariable('no_cont_ub', 'f4', ( 'sector','time','lat','lon'))\n",
    "value[:,:,:,:]=no_cont_ub\n",
    "# lats[:]=latgc\n",
    "# lons[:]=longc\n",
    "ds.close()\n",
    "# pm_cont = np.transpose(Dataset(out_dir +'LOS_PM_NO_O3.nc', mode='r').variables['pm_cont'][:])\n",
    "# pm_cont_lb = np.transpose(Dataset(out_dir +'LOS_PM_NO_O3.nc', mode='r').variables['pm_cont_lb'][:])\n",
    "# pm_cont_ub = np.transpose(Dataset(out_dir +'LOS_PM_NO_O3.nc', mode='r').variables['pm_cont_ub'][:])\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "        \n",
    "# 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",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "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",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "# 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": "35bcb2c9",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}