{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import pickle, pandas as pd, pmdarima as pmd, geopandas as gpd,numpy as np\n", "from random import choice\n", "import crf_utilities_again as crf_util\n", "import xarray as xr, geopandas as gpd, rioxarray, os\n", "from shapely.geometry import box, mapping\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "regions_all = gpd.read_file(\n", " '/home/gchossie/2020_aq/monitor_based/data/all_regions.geojson')" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "results_geo = regions_all.copy()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results_geo" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "for col in ['NO2 mean', 'NO2 2.5', 'NO2 97.5',\n", " 'PM2.5 mean', 'PM2.5 2.5', 'PM2.5 97.5',\n", " 'O3 mean', 'O3 2.5', 'O3 97.5']:\n", " results_geo.loc[:,col] = 0" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "results_geo.reset_index(drop=True,inplace=True)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
EnsembleCountryRegionPopulationgeometryNO2 meanNO2 2.5NO2 97.5PM2.5 meanPM2.5 2.5PM2.5 97.5O3 meanO3 2.5O3 97.5
0JapanJapanAichi7220143.5MULTIPOLYGON (((137.09523 34.65330, 137.09546 ...000000000
1JapanJapanOsaka8177242.0MULTIPOLYGON (((135.31573 34.43310, 135.31544 ...000000000
2JapanJapanShizuoka3460957.0MULTIPOLYGON (((141.46085 24.24579, 141.46115 ...000000000
\n", "
" ], "text/plain": [ " Ensemble Country Region Population \\\n", "0 Japan Japan Aichi 7220143.5 \n", "1 Japan Japan Osaka 8177242.0 \n", "2 Japan Japan Shizuoka 3460957.0 \n", "\n", " geometry NO2 mean NO2 2.5 \\\n", "0 MULTIPOLYGON (((137.09523 34.65330, 137.09546 ... 0 0 \n", "1 MULTIPOLYGON (((135.31573 34.43310, 135.31544 ... 0 0 \n", "2 MULTIPOLYGON (((141.46085 24.24579, 141.46115 ... 0 0 \n", "\n", " NO2 97.5 PM2.5 mean PM2.5 2.5 PM2.5 97.5 O3 mean O3 2.5 O3 97.5 \n", "0 0 0 0 0 0 0 0 \n", "1 0 0 0 0 0 0 0 \n", "2 0 0 0 0 0 0 0 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results_geo.head(3)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "raw_dist = {}" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "N = 10000\n", "samples = np.random.normal(size=N)\n", "for idx in results_geo.index:\n", " country = results_geo.loc[idx,'Country']\n", " reg = results_geo.loc[idx,'Region']\n", " if country not in raw_dist.keys():\n", " raw_dist[country] = {reg: {}}\n", " else:\n", " raw_dist[country] = {**raw_dist[country], **{reg: {}}}\n", " for sp in ['O3', 'PM2.5', 'NO2']:\n", " raw_dist[country][reg][sp] = {}\n", " for cause in ['Lower respiratory infections',\n", " 'Upper respiratory infections',\n", " 'Chronic respiratory diseases',\n", " 'Cardiovascular diseases',\n", " 'All causes']:\n", " mort = crf_util.get_mort(country, reg, sp, cause='All causes', scen='', \n", " samples=samples, method='kriging', sens=True, \n", " counterfactual=False)\n", " if mort is None: \n", " raw_dist[country][reg][sp][cause] = 0 \n", " continue\n", " else:\n", " raw_dist[country][reg][sp][cause] = mort\n", " try:\n", " results_geo.loc[idx,f'{sp} mean'] = mort.mean()\n", " results_geo.loc[idx,f'{sp} 2.5'] = np.percentile(mort,2.5)\n", " results_geo.loc[idx,f'{sp} 97.5'] = np.percentile(mort,97.5)\n", " except: pass" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results_geo.loc[results_geo.Region == 'Massachusetts',:]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results_geo['NO2 mean'].sum() + results_geo['PM2.5 mean'].sum() + results_geo['O3 mean'].sum() " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ctr = 'United States'; sp = 'PM2.5'\n", "regs = list(set(results_geo.loc[results_geo.Country == ctr, 'Region']))\n", "np.percentile(np.sum([raw_dist[ctr][reg][sp] for reg in regs], axis=0),[2.5, 50, 97.5])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results_geo.loc[results_geo.Country == ctr, 'NO2 mean'].sum()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sp = 'PM2.5'\n", "np.percentile(np.sum([raw_dist[results_geo.loc[idx,'Country']][results_geo.loc[idx,'Region']][sp] for idx in results_geo.index], axis=0),50)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "results_geo.to_file('./results/result_moni_mort_oct_15.geojson',driver='GeoJSON')\n", "#results_geo.to_file('./results/result_moni_mort_jul9_counterfactual.geojson',driver='GeoJSON')" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "with open('./results/raw_mort_dist_oct_21_sens.pkl','wb') as f:\n", " pickle.dump(raw_dist,f)\n", "#with open('./results/raw_mort_dist_jul9_counterfactual.pkl','wb') as f:\n", "# pickle.dump(raw_dist,f)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Freight US", "language": "python", "name": "freight_us" }, "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.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }