diff --git a/analysis/a00_seasonal_forecast_analysis.ipynb b/analysis/a00_seasonal_forecast_analysis.ipynb index b0b0e7f..f277f19 100644 --- a/analysis/a00_seasonal_forecast_analysis.ipynb +++ b/analysis/a00_seasonal_forecast_analysis.ipynb @@ -39,20 +39,18 @@ ], "source": [ "import os\n", - "from pathlib import Path\n", "from dotenv import load_dotenv\n", "import xarray as xr\n", "\n", "import geopandas as gpd\n", - "import pandas as pd\n", "import numpy as np\n", - "import math\n", "import matplotlib.pyplot as plt\n", "\n", "import analysis_utils\n", "import isku_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -78,7 +76,7 @@ "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", "\n", - "#Climate Data\n", + "# Climate Data\n", "TAS_FORECAST_URI = os.environ[\"POREALLAS_TAS_FORECAST_URI\"]\n", "ERA5_URI = os.environ[\"POREALLAS_ERA5_URI\"]\n", "TAS_FORECAST_RAW = os.environ[\"POREALLAS_TAS_FORECAST_RAW_URI\"]" @@ -91,10 +89,10 @@ "metadata": {}, "outputs": [], "source": [ - "#Projection Effects\n", - "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated = False)\n", - "baseline_period = analysis_utils.get_baseline_period(effect, years = 30)\n", - "#Impact Regions\n", + "# Projection Effects\n", + "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated=False)\n", + "baseline_period = analysis_utils.get_baseline_period(effect, years=30)\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(os.path.join(DATA_DIR, IMPACT_REGION_POLYGONS))\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -102,9 +100,11 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(os.path.join(DATA_DIR, SOCIOECONOMICS_URI))\n", - "socioeconomics = socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop', 'gdppc', 'iso3']]" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(os.path.join(DATA_DIR, SOCIOECONOMICS_URI))\n", + "socioeconomics = socioeconomics.sel(year=2026)[\n", + " [\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\", \"gdppc\", \"iso3\"]\n", + "]" ] }, { @@ -122,11 +122,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Compute the Mortality Impact relative to a Baseline Period for the Projected Effects\n", - "impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = True, baseline_period=baseline_period, hotonly = False, rate = True, age_weight = True)\n", + "# Compute the Mortality Impact relative to a Baseline Period for the Projected Effects\n", + "impact = analysis_utils.compute_impact(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=False,\n", + " rate=True,\n", + " age_weight=True,\n", + ")\n", "\n", "##For a single age-cohort, use:\n", - "#impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = False, baseline_period=baseline_period, hotonly = False, rate = True, age_weight = False, cohort = 'age65plus')" + "# impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = False, baseline_period=baseline_period, hotonly = False, rate = True, age_weight = False, cohort = 'age65plus')" ] }, { @@ -169,7 +177,9 @@ } ], "source": [ - "_polygons_impact = analysis_utils.xarray_to_gpd(impact.mean(dim = ['number', 'sample']), _polygons)" + "_polygons_impact = analysis_utils.xarray_to_gpd(\n", + " impact.mean(dim=[\"number\", \"sample\"]), _polygons\n", + ")" ] }, { @@ -193,12 +203,12 @@ "# Plot Impact on a Monthly Basis\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " sup_title=f\"\",\n", - " cbar_label= impact.units,\n", - " month_order= [8, 9, 10, 11, 12, 1]\n", - ")\n" + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " sup_title=\"\",\n", + " cbar_label=impact.units,\n", + " month_order=[8, 9, 10, 11, 12, 1],\n", + ")" ] }, { @@ -222,11 +232,11 @@ "# Plot 6-month Cumulative Impact\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " agg = 'sum',\n", - " cm = 'bwr',\n", - " sup_title=f\"\",\n", - " cbar_label= impact.units,\n", + " col=\"age_weighted_impact\",\n", + " agg=\"sum\",\n", + " cm=\"bwr\",\n", + " sup_title=\"\",\n", + " cbar_label=impact.units,\n", ")" ] }, @@ -248,14 +258,14 @@ } ], "source": [ - "#Ex: Plot 6-month SD\n", + "# Ex: Plot 6-month SD\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " agg = 'std',\n", - " cm = 'YlOrRd',\n", - " sup_title=f\"\",\n", - " cbar_label= impact.units,\n", + " col=\"age_weighted_impact\",\n", + " agg=\"std\",\n", + " cm=\"YlOrRd\",\n", + " sup_title=\"\",\n", + " cbar_label=impact.units,\n", ")" ] }, @@ -298,7 +308,7 @@ "source": [ "forecast = xr.open_zarr(os.path.join(DATA_DIR, TAS_FORECAST_URI))\n", "# Gridded to Impact Region\n", - "forecast_ir = isku_utils.grid_to_ir(forecast.mean(dim = 'number'), savefile = None)" + "forecast_ir = isku_utils.grid_to_ir(forecast.mean(dim=\"number\"), savefile=None)" ] }, { @@ -308,7 +318,9 @@ "metadata": {}, "outputs": [], "source": [ - "_polygons_forecast = analysis_utils.xarray_to_gpd(forecast_ir['tas'].groupby('time.month').mean(), _polygons)" + "_polygons_forecast = analysis_utils.xarray_to_gpd(\n", + " forecast_ir[\"tas\"].groupby(\"time.month\").mean(), _polygons\n", + ")" ] }, { @@ -331,13 +343,13 @@ "source": [ "fig = analysis_utils.plot_monthly(\n", " _polygons_forecast,\n", - " col='tas',\n", - " vmin = 273.15,\n", - " vmax = 303,\n", - " cm = 'Reds',\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS [K]',\n", - " month_order= [8, 9, 10, 11, 12, 1]\n", + " col=\"tas\",\n", + " vmin=273.15,\n", + " vmax=303,\n", + " cm=\"Reds\",\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS [K]\",\n", + " month_order=[8, 9, 10, 11, 12, 1],\n", ")" ] }, @@ -356,10 +368,15 @@ "metadata": {}, "outputs": [], "source": [ - "#Raw Forecast/Hindcast\n", + "# Raw Forecast/Hindcast\n", "s51 = xr.open_zarr(os.path.join(DATA_DIR, TAS_FORECAST_RAW))\n", - "forecast = s51.sel(time=slice('2026-08-01', '2027-01-31')).groupby(\"time.month\").mean()\n", - "hindcast = s51.sel(time=slice('1993-08-01', '2016-01-31')).groupby(\"time.month\").mean().sel(month = [1, 8, 9, 10, 11, 12])" + "forecast = s51.sel(time=slice(\"2026-08-01\", \"2027-01-31\")).groupby(\"time.month\").mean()\n", + "hindcast = (\n", + " s51.sel(time=slice(\"1993-08-01\", \"2016-01-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[1, 8, 9, 10, 11, 12])\n", + ")" ] }, { @@ -369,7 +386,7 @@ "metadata": {}, "outputs": [], "source": [ - "anomaly = forecast['tas']-hindcast['tas']\n", + "anomaly = forecast[\"tas\"] - hindcast[\"tas\"]\n", "anomaly.name = \"forecast_anomaly\"" ] }, @@ -394,8 +411,10 @@ } ], "source": [ - "anomaly_ir = isku_utils.grid_to_ir(anomaly.compute(), savefile = None)\n", - "_polygons_anomaly = analysis_utils.xarray_to_gpd(anomaly_ir.mean(dim = 'number'), _polygons)" + "anomaly_ir = isku_utils.grid_to_ir(anomaly.compute(), savefile=None)\n", + "_polygons_anomaly = analysis_utils.xarray_to_gpd(\n", + " anomaly_ir.mean(dim=\"number\"), _polygons\n", + ")" ] }, { @@ -418,11 +437,11 @@ "source": [ "fig = analysis_utils.plot_monthly(\n", " _polygons_anomaly,\n", - " col='value',\n", - " cm = 'bwr',\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS Anomaly [C]',\n", - " month_order= [8, 9, 10, 11, 12, 1]\n", + " col=\"value\",\n", + " cm=\"bwr\",\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS Anomaly [C]\",\n", + " month_order=[8, 9, 10, 11, 12, 1],\n", ")" ] }, @@ -444,14 +463,14 @@ } ], "source": [ - "#6 month average anomaly\n", + "# 6 month average anomaly\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_anomaly,\n", - " col='value',\n", - " agg = 'mean',\n", - " cm = 'bwr',\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS Anomaly [C]',\n", + " col=\"value\",\n", + " agg=\"mean\",\n", + " cm=\"bwr\",\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS Anomaly [C]\",\n", ")" ] }, @@ -470,9 +489,15 @@ "metadata": {}, "outputs": [], "source": [ - "analysis_utils.make_csv(effect, socioeconomics, ensemble = True, \n", - " baseline_period=baseline_period, hotonly = \"coldonly\", \n", - " rate = False, age_weight = True)" + "analysis_utils.make_csv(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=\"coldonly\",\n", + " rate=False,\n", + " age_weight=True,\n", + ")" ] }, { @@ -492,10 +517,14 @@ "source": [ "### TEST DATA ###\n", "# IR Climate Data\n", - "forecast_gmfd = xr.open_zarr(os.path.join(DATA_DIR, \"daily_forecast_GMFD_tas.zarr\"))-273.15\n", - "reanalysis_gmfd = xr.open_zarr(os.path.join(DATA_DIR, \"daily_reanalysis_GMFD_tas.zarr\"))-273.15\n", + "forecast_gmfd = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"daily_forecast_GMFD_tas.zarr\")) - 273.15\n", + ")\n", + "reanalysis_gmfd = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"daily_reanalysis_GMFD_tas.zarr\")) - 273.15\n", + ")\n", "\n", - "#Intermediate Betas/MMT from Project.Py\n", + "# Intermediate Betas/MMT from Project.Py\n", "betas_mmt = xr.open_zarr(os.path.join(DATA_DIR, \"beta_mmt_GMFD.zarr\"))" ] }, @@ -506,21 +535,15 @@ "metadata": {}, "outputs": [], "source": [ - "import matplotlib.patheffects as pe\n", - "from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec\n", - "\n", - "import cartopy.crs as ccrs\n", - "from scipy.interpolate import interp1d\n", - "from shapely import box\n", - "\n", "import diagnostic_utils\n", + "\n", "importlib.reload(diagnostic_utils)\n", "\n", "### Analysis Options ###\n", - "rate = True # Mortality Rate vs Total Deaths\n", - "hotonly = False # Hot-Only Deaths\n", - "age_weight = False #Proportional Age-weighting vs single age group\n", - "density_hist = False #Controls daily data histograms (not computation normalization)\n", + "rate = True # Mortality Rate vs Total Deaths\n", + "hotonly = False # Hot-Only Deaths\n", + "age_weight = False # Proportional Age-weighting vs single age group\n", + "density_hist = False # Controls daily data histograms (not computation normalization)\n", "\n", "### Region Selection ###\n", "region_filters = [\n", @@ -542,27 +565,38 @@ "\n", "# Compute Impact\n", "projected = effect\n", - "regional_sum = analysis_utils.compute_impact(projected, socioeconomics, baseline_period= baseline_period, hotonly = hotonly, rate = rate, age_weight = age_weight)\n", + "regional_sum = analysis_utils.compute_impact(\n", + " projected,\n", + " socioeconomics,\n", + " baseline_period=baseline_period,\n", + " hotonly=hotonly,\n", + " rate=rate,\n", + " age_weight=age_weight,\n", + ")\n", "\n", - "#Get Months\n", + "# Get Months\n", "months = regional_sum.month.values\n", "\n", "### Load Climate Data ###\n", - "#Select representative ensemble member\n", - "forecast_slice = forecast_gmfd['tas'].sel(number = 1)\n", - "forecast_local = forecast_slice.sel(region = region_filter)\n", + "# Select representative ensemble member\n", + "forecast_slice = forecast_gmfd[\"tas\"].sel(number=1)\n", + "forecast_local = forecast_slice.sel(region=region_filter)\n", "\n", - "#Select representative year\n", - "reanalysis_slice = reanalysis_gmfd['tas'].sel(time=slice('2000-01-01', '2000-12-31'))\n", + "# Select representative year\n", + "reanalysis_slice = reanalysis_gmfd[\"tas\"].sel(time=slice(\"2000-01-01\", \"2000-12-31\"))\n", "month_mask = reanalysis_slice.month.isin(months).compute()\n", "reanalysis_slice = reanalysis_slice.where(month_mask, drop=True)\n", - "reanalysis_local = reanalysis_slice.sel(region = region_filter)\n", + "reanalysis_local = reanalysis_slice.sel(region=region_filter)\n", "\n", - "combined = np.concat([forecast_local.values.flatten(), reanalysis_local.values.flatten()])\n", + "combined = np.concat(\n", + " [forecast_local.values.flatten(), reanalysis_local.values.flatten()]\n", + ")\n", "xmin, xmax = np.percentile(combined, [1, 100])\n", "\n", "# Compute Cumulative Effect\n", - "cdf_data, max_cdf = diagnostic_utils.compute_cumulative_effect(forecast_slice, reanalysis_slice, region_filter, months, hotonly = hotonly)" + "cdf_data, max_cdf = diagnostic_utils.compute_cumulative_effect(\n", + " forecast_slice, reanalysis_slice, region_filter, months, hotonly=hotonly\n", + ")" ] }, { @@ -583,35 +617,63 @@ } ], "source": [ - "#Histogram\n", + "# Histogram\n", "\n", "b_color = \"#999999\"\n", "f_color = \"#ffbb6f\"\n", "r_color = \"#5e4c5f\"\n", "\n", "fig, (ax_hist, ax_cdf) = plt.subplots(\n", - " 2, 1, sharex=True, figsize=(6, 6),\n", - " gridspec_kw={'height_ratios': [3, 1], 'hspace': 0.08}\n", + " 2,\n", + " 1,\n", + " sharex=True,\n", + " figsize=(6, 6),\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.08},\n", ")\n", "\n", "\n", - "ref_vals = betas_mmt['mmt'].sel(region=region_filter).sel(age_cohort='age65plus')\n", - "betas = betas_mmt['beta_hotonly'] if hotonly else betas_mmt['beta']\n", - "da_temp_bins = betas.sel(region=region_filter).sel(age_cohort='age65plus')\n", - "da_temp_bins['tas_bin'] = da_temp_bins['tas_bin'] \n", + "ref_vals = betas_mmt[\"mmt\"].sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + "betas = betas_mmt[\"beta_hotonly\"] if hotonly else betas_mmt[\"beta\"]\n", + "da_temp_bins = betas.sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + "da_temp_bins[\"tas_bin\"] = da_temp_bins[\"tas_bin\"]\n", "\n", "# Density histograms on hidden twin axis\n", - "bins = da_temp_bins.coords['tas_bin']\n", - "ax_hist.hist(reanalysis_local.values.flatten(), bins=bins, alpha=0.5, label=\"ERA5 (2000)\", density=density_hist, color=b_color)\n", - "ax_hist.hist(forecast_local.values.flatten(), bins=bins, alpha=0.5, label=\"Forecast\", density=density_hist, color=f_color)\n", - "ax_hist.axvline(forecast_local.values.flatten().mean(), color=f_color, linewidth=1, linestyle='--')\n", - "ax_hist.axvline(reanalysis_local.values.flatten().mean(), color=b_color, linewidth=1, linestyle='--')\n", + "bins = da_temp_bins.coords[\"tas_bin\"]\n", + "ax_hist.hist(\n", + " reanalysis_local.values.flatten(),\n", + " bins=bins,\n", + " alpha=0.5,\n", + " label=\"ERA5 (2000)\",\n", + " density=density_hist,\n", + " color=b_color,\n", + ")\n", + "ax_hist.hist(\n", + " forecast_local.values.flatten(),\n", + " bins=bins,\n", + " alpha=0.5,\n", + " label=\"Forecast\",\n", + " density=density_hist,\n", + " color=f_color,\n", + ")\n", + "ax_hist.axvline(\n", + " forecast_local.values.flatten().mean(), color=f_color, linewidth=1, linestyle=\"--\"\n", + ")\n", + "ax_hist.axvline(\n", + " reanalysis_local.values.flatten().mean(), color=b_color, linewidth=1, linestyle=\"--\"\n", + ")\n", "ax_hist.set_ylabel(\"Count [Days]\")\n", "ax_hist.set_xlim(xmin, xmax)\n", "\n", "# Response curve on the primary (left) axis\n", "ax_hist2 = ax_hist.twinx()\n", - "ax_hist2.step(da_temp_bins.coords['tas_bin'], da_temp_bins.values, where='pre', color=r_color, linewidth=2, label=\"Response curve\")\n", + "ax_hist2.step(\n", + " da_temp_bins.coords[\"tas_bin\"],\n", + " da_temp_bins.values,\n", + " where=\"pre\",\n", + " color=r_color,\n", + " linewidth=2,\n", + " label=\"Response curve\",\n", + ")\n", "ax_hist2.axvline(ref_vals.values, color=r_color, linewidth=1, label=\"MMT\")\n", "ax_hist2.set_ylabel(\"Daily Response 65+ \\n [deaths/100k]\")\n", "ax_hist2.set_ylim(0, 10)\n", @@ -621,23 +683,36 @@ "era5_centers, era5_cdf, fc_centers, fc_cdf = cdf_data\n", "\n", "ax_cdf.set_ylabel(\"6-Month Cum.\\n [Deaths/100k]\")\n", - "ax_cdf.set_ylim(0, max_cdf+2)\n", + "ax_cdf.set_ylim(0, max_cdf + 2)\n", "ax_cdf.set_xlim(xmin, xmax)\n", "ax_cdf.set_xlabel(\"Tas [C]\")\n", - "ax_cdf.step(fc_centers, fc_cdf, where='post', color=f_color, label=\"Forecast CDF\")\n", - "ax_cdf.step(era5_centers, era5_cdf, where='post', color=b_color, label=\"ERA5 CDF\")\n", + "ax_cdf.step(fc_centers, fc_cdf, where=\"post\", color=f_color, label=\"Forecast CDF\")\n", + "ax_cdf.step(era5_centers, era5_cdf, where=\"post\", color=b_color, label=\"ERA5 CDF\")\n", "\n", - "#Legend\n", - "handles = (ax_hist.get_legend_handles_labels()[0]\n", - " + ax_hist2.get_legend_handles_labels()[0]\n", - " + ax_cdf.get_legend_handles_labels()[0])\n", - "labels = (ax_hist.get_legend_handles_labels()[1]\n", - " + ax_hist2.get_legend_handles_labels()[1]\n", - " + ax_cdf.get_legend_handles_labels()[1])\n", - "ax_cdf.legend(handles, labels, fontsize=7, loc='upper center',\n", - " bbox_to_anchor=(0.5, -0.35), ncol=3, frameon=False)\n", - "ax_hist.set_title(f'{region_filter}')\n", - "fig.savefig(f'response_diagnostic_{region_filter}_2607.png', dpi=600, bbox_inches=\"tight\")" + "# Legend\n", + "handles = (\n", + " ax_hist.get_legend_handles_labels()[0]\n", + " + ax_hist2.get_legend_handles_labels()[0]\n", + " + ax_cdf.get_legend_handles_labels()[0]\n", + ")\n", + "labels = (\n", + " ax_hist.get_legend_handles_labels()[1]\n", + " + ax_hist2.get_legend_handles_labels()[1]\n", + " + ax_cdf.get_legend_handles_labels()[1]\n", + ")\n", + "ax_cdf.legend(\n", + " handles,\n", + " labels,\n", + " fontsize=7,\n", + " loc=\"upper center\",\n", + " bbox_to_anchor=(0.5, -0.35),\n", + " ncol=3,\n", + " frameon=False,\n", + ")\n", + "ax_hist.set_title(f\"{region_filter}\")\n", + "fig.savefig(\n", + " f\"response_diagnostic_{region_filter}_2607.png\", dpi=600, bbox_inches=\"tight\"\n", + ")" ] }, { diff --git a/analysis/a01_august_seasonal_forecast_results.ipynb b/analysis/a01_august_seasonal_forecast_results.ipynb index c26099c..3f9ed01 100644 --- a/analysis/a01_august_seasonal_forecast_results.ipynb +++ b/analysis/a01_august_seasonal_forecast_results.ipynb @@ -38,20 +38,19 @@ ], "source": [ "import os\n", - "from pathlib import Path\n", - "from dotenv import load_dotenv, find_dotenv\n", + "from dotenv import load_dotenv\n", "import xarray as xr\n", "\n", "import geopandas as gpd\n", "import pandas as pd\n", "import numpy as np\n", - "import math\n", "import matplotlib.pyplot as plt\n", "\n", "import analysis_utils\n", "import isku_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -63,7 +62,7 @@ "metadata": {}, "outputs": [], "source": [ - "adj = xr.open_zarr('/home/emily_zuetell/projects/poreallas/data/m1_forecast_adj.zarr')" + "adj = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/m1_forecast_adj.zarr\")" ] }, { @@ -73,7 +72,9 @@ "metadata": {}, "outputs": [], "source": [ - "old = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/forecast_adj.zarr\").compute()" + "old = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/forecast_adj.zarr\"\n", + ").compute()" ] }, { @@ -115,13 +116,15 @@ ], "source": [ "fig, ax = plt.subplots()\n", - "old_qdm = old['tas'].sortby(\"time\").sel(time = slice(\"2027-01-01\", \"2027-01-31\"))\n", - "old_qdm = old_qdm.where(old_qdm<1000, drop = True)\n", - "new_qdm = adj['tas']\n", + "old_qdm = old[\"tas\"].sortby(\"time\").sel(time=slice(\"2027-01-01\", \"2027-01-31\"))\n", + "old_qdm = old_qdm.where(old_qdm < 1000, drop=True)\n", + "new_qdm = adj[\"tas\"]\n", "\n", - "sns.histplot(old_qdm.values.flatten(), kde=True, bins=30, ax = ax, stat = 'density')\n", - "sns.histplot(new_qdm.values.flatten(), kde=True, bins=30, ax = ax, stat = 'density', color = 'k')\n", - "ax.set_xlabel('Adj_Forecast [K]')" + "sns.histplot(old_qdm.values.flatten(), kde=True, bins=30, ax=ax, stat=\"density\")\n", + "sns.histplot(\n", + " new_qdm.values.flatten(), kde=True, bins=30, ax=ax, stat=\"density\", color=\"k\"\n", + ")\n", + "ax.set_xlabel(\"Adj_Forecast [K]\")" ] }, { @@ -141,12 +144,14 @@ "source": [ "load_dotenv()\n", "DATA_DIR = os.environ[\"DATA_DIR\"]\n", - "#EFFECTS_URI = os.environ[\"POREALLAS_EFFECTS_URI\"]\n", - "EFFECTS_URI = \"/home/emily_zuetell/projects/poreallas/data/2608_effect_gamma_sampled.zarr\"\n", + "# EFFECTS_URI = os.environ[\"POREALLAS_EFFECTS_URI\"]\n", + "EFFECTS_URI = (\n", + " \"/home/emily_zuetell/projects/poreallas/data/2608_effect_gamma_sampled.zarr\"\n", + ")\n", "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", "\n", - "#Climate Data\n", + "# Climate Data\n", "TAS_FORECAST_URI = os.environ[\"POREALLAS_TAS_FORECAST_URI\"]\n", "ERA5_URI = os.environ[\"POREALLAS_ERA5_URI\"]" ] @@ -158,10 +163,10 @@ "metadata": {}, "outputs": [], "source": [ - "#Projection Effects\n", - "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated = False)\n", - "baseline_period = analysis_utils.get_baseline_period(effect, years = 30)\n", - "#Impact Regions\n", + "# Projection Effects\n", + "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated=False)\n", + "baseline_period = analysis_utils.get_baseline_period(effect, years=30)\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(IMPACT_REGION_POLYGONS)\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -169,9 +174,11 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", - "socioeconomics = socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop', 'gdppc', 'iso3']]" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", + "socioeconomics = socioeconomics.sel(year=2026)[\n", + " [\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\", \"gdppc\", \"iso3\"]\n", + "]" ] }, { @@ -189,11 +196,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Compute the Mortality Impact relative to a Baseline Period for the Projected Effects\n", - "impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = True, baseline_period=baseline_period, hotonly = \"hotonly\", rate = False, age_weight = True)\n", + "# Compute the Mortality Impact relative to a Baseline Period for the Projected Effects\n", + "impact = analysis_utils.compute_impact(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=\"hotonly\",\n", + " rate=False,\n", + " age_weight=True,\n", + ")\n", "\n", "##For a single age-cohort, use:\n", - "#impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = False, baseline_period=baseline_period, hotonly = False, rate = True, age_weight = False, cohort = 'age65plus')" + "# impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = False, baseline_period=baseline_period, hotonly = False, rate = True, age_weight = False, cohort = 'age65plus')" ] }, { @@ -203,7 +218,7 @@ "metadata": {}, "outputs": [], "source": [ - "#Total by Country\n", + "# Total by Country\n", "impact = impact.assign_coords(iso3=socioeconomics[\"iso3\"]).compute()" ] }, @@ -214,7 +229,7 @@ "metadata": {}, "outputs": [], "source": [ - "impact_country = impact.groupby('iso3').sum()" + "impact_country = impact.groupby(\"iso3\").sum()" ] }, { @@ -225,16 +240,24 @@ "outputs": [], "source": [ "rate = False\n", - "hotonly = 'hotonly'\n", + "hotonly = \"hotonly\"\n", "\n", "rate_l = \"rate\" if rate else \"total\"\n", "\n", "ds = impact_country\n", "# Monthly Stats\n", - "stat_cols = ['median', 'p17', 'p83', 'likely_range_IPCC', \n", - " 'mean','std', \n", - " 'min', 'max', \n", - " 'p10', 'p90']\n", + "stat_cols = [\n", + " \"median\",\n", + " \"p17\",\n", + " \"p83\",\n", + " \"likely_range_IPCC\",\n", + " \"mean\",\n", + " \"std\",\n", + " \"min\",\n", + " \"max\",\n", + " \"p10\",\n", + " \"p90\",\n", + "]\n", "stats_impact = analysis_utils.compute_stats(ds, dim=[\"number\", \"sample\"])" ] }, @@ -246,7 +269,7 @@ "outputs": [], "source": [ "stats_df = stats_impact.to_dataframe().reset_index()\n", - "wide = stats_df.pivot(index=['iso3'], columns=\"month\", values=stat_cols)" + "wide = stats_df.pivot(index=[\"iso3\"], columns=\"month\", values=stat_cols)" ] }, { @@ -258,8 +281,7 @@ "source": [ "wide.columns = [f\"month {m} {stat}\" for stat, m in wide.columns]\n", "wide = wide.reset_index()\n", - "wide.to_csv(f\"2608_{hotonly}_{rate_l}_all_stats_bycountry.csv\", index = False)\n", - "\n" + "wide.to_csv(f\"2608_{hotonly}_{rate_l}_all_stats_bycountry.csv\", index=False)" ] }, { @@ -291,15 +313,25 @@ "outputs": [], "source": [ "# 6-month stats\n", - "mo6 = ds.sum(dim = 'month')\n", + "mo6 = ds.sum(dim=\"month\")\n", "stats_mo6 = analysis_utils.compute_stats(mo6, dim=[\"number\", \"sample\"])\n", "mo6_df = stats_mo6.to_dataframe().reset_index()\n", - "mo6_out = mo6_df[[\"iso3\", \n", - " 'median', 'p17', 'p83', 'likely_range_IPCC', \n", - " 'mean','std', \n", - " 'min', 'max', \n", - " 'p10', 'p90']]\n", - "mo6_out.to_csv(f\"2608_{hotonly}_6mo_{rate_l}_all_stats_bycountry.csv\", index = False)" + "mo6_out = mo6_df[\n", + " [\n", + " \"iso3\",\n", + " \"median\",\n", + " \"p17\",\n", + " \"p83\",\n", + " \"likely_range_IPCC\",\n", + " \"mean\",\n", + " \"std\",\n", + " \"min\",\n", + " \"max\",\n", + " \"p10\",\n", + " \"p90\",\n", + " ]\n", + "]\n", + "mo6_out.to_csv(f\"2608_{hotonly}_6mo_{rate_l}_all_stats_bycountry.csv\", index=False)" ] }, { @@ -320,9 +352,13 @@ "outputs": [], "source": [ "# Monthly Global Average (Net)\n", - "global_avg_month = impact.sum(dim = ['region']).quantile([0.17, 0.5, 0.83], dim = ['number', 'sample']).to_dataframe()\n", + "global_avg_month = (\n", + " impact.sum(dim=[\"region\"])\n", + " .quantile([0.17, 0.5, 0.83], dim=[\"number\", \"sample\"])\n", + " .to_dataframe()\n", + ")\n", "table = global_avg_month.unstack(level=\"quantile\")\n", - "table.to_csv('global_avg_month_net.csv')" + "table.to_csv(\"global_avg_month_net.csv\")" ] }, { @@ -333,9 +369,13 @@ "outputs": [], "source": [ "# Monthly Global Average (Hot Only)\n", - "global_avg_month = impact.sum(dim = ['region']).quantile([0.17, 0.5, 0.83], dim = ['number', 'sample']).to_dataframe()\n", + "global_avg_month = (\n", + " impact.sum(dim=[\"region\"])\n", + " .quantile([0.17, 0.5, 0.83], dim=[\"number\", \"sample\"])\n", + " .to_dataframe()\n", + ")\n", "table = global_avg_month.unstack(level=\"quantile\")\n", - "table.to_csv('global_avg_month_hotonly.csv')" + "table.to_csv(\"global_avg_month_hotonly.csv\")" ] }, { @@ -354,7 +394,7 @@ } ], "source": [ - "_polygons_impact = analysis_utils.xarray_to_gpd(impact.mean(dim = 'number'), _polygons)" + "_polygons_impact = analysis_utils.xarray_to_gpd(impact.mean(dim=\"number\"), _polygons)" ] }, { @@ -378,10 +418,10 @@ "# Plot Impact on a Monthly Basis\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " sup_title=f\"\",\n", - " cbar_label= impact.units,\n", + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " sup_title=\"\",\n", + " cbar_label=impact.units,\n", ")" ] }, @@ -395,11 +435,11 @@ "# Plot 6-month Cumulative Impact\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " agg = 'sum',\n", - " cm = 'bwr',\n", - " sup_title=f\"\",\n", - " cbar_label= impact.units,\n", + " col=\"age_weighted_impact\",\n", + " agg=\"sum\",\n", + " cm=\"bwr\",\n", + " sup_title=\"\",\n", + " cbar_label=impact.units,\n", ")" ] }, @@ -410,14 +450,14 @@ "metadata": {}, "outputs": [], "source": [ - "#Ex: Plot 6-month SD\n", + "# Ex: Plot 6-month SD\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " agg = 'std',\n", - " cm = 'YlOrRd',\n", - " sup_title=f\"\",\n", - " cbar_label= impact.units,\n", + " col=\"age_weighted_impact\",\n", + " agg=\"std\",\n", + " cm=\"YlOrRd\",\n", + " sup_title=\"\",\n", + " cbar_label=impact.units,\n", ")" ] }, @@ -446,7 +486,7 @@ "source": [ "forecast = xr.open_zarr(TAS_FORECAST_URI)\n", "# Gridded to Impact Region\n", - "forecast_ir = isku_utils.grid_to_ir(forecast.mean(dim = 'number'), savefile = None)" + "forecast_ir = isku_utils.grid_to_ir(forecast.mean(dim=\"number\"), savefile=None)" ] }, { @@ -456,7 +496,7 @@ "metadata": {}, "outputs": [], "source": [ - "_polygons_forecast = analysis_utils.xarray_to_gpd(forecast_ir['tas'], _polygons)" + "_polygons_forecast = analysis_utils.xarray_to_gpd(forecast_ir[\"tas\"], _polygons)" ] }, { @@ -468,12 +508,12 @@ "source": [ "fig = analysis_utils.plot_monthly(\n", " _polygons_forecast,\n", - " col='tas',\n", - " vmin = 273.15,\n", - " vmax = 303,\n", - " cm = 'Reds',\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS [K]',\n", + " col=\"tas\",\n", + " vmin=273.15,\n", + " vmax=303,\n", + " cm=\"Reds\",\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS [K]\",\n", ")" ] }, @@ -492,10 +532,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Raw Forecast/Hindcast\n", - "s51 = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/parsed/08_forecast_parsed.zarr\")\n", - "forecast_raw = s51.sel(time=slice('2026-08-01', '2027-01-31')).groupby(\"time.month\").mean()\n", - "hindcast_raw = s51.sel(time=slice('1993-08-01', '2016-01-31')).groupby(\"time.month\").mean().sel(month = [1, 8, 9, 10, 11, 12])" + "# Raw Forecast/Hindcast\n", + "s51 = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/parsed/08_forecast_parsed.zarr\"\n", + ")\n", + "forecast_raw = (\n", + " s51.sel(time=slice(\"2026-08-01\", \"2027-01-31\")).groupby(\"time.month\").mean()\n", + ")\n", + "hindcast_raw = (\n", + " s51.sel(time=slice(\"1993-08-01\", \"2016-01-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[1, 8, 9, 10, 11, 12])\n", + ")" ] }, { @@ -505,12 +554,12 @@ "metadata": {}, "outputs": [], "source": [ - "#Compute monthly average anomaly\n", - "anomaly = forecast_raw['tas']-hindcast_raw['tas']\n", + "# Compute monthly average anomaly\n", + "anomaly = forecast_raw[\"tas\"] - hindcast_raw[\"tas\"]\n", "anomaly.name = \"forecast_anomaly\"\n", "\n", "## Gridded to Impact Region\n", - "anomaly_ir = isku_utils.grid_to_ir(anomaly.mean(dim = 'number'), savefile = None)\n", + "anomaly_ir = isku_utils.grid_to_ir(anomaly.mean(dim=\"number\"), savefile=None)\n", "_polygons_anomaly = analysis_utils.xarray_to_gpd(anomaly_ir, _polygons)" ] }, @@ -523,14 +572,15 @@ "source": [ "fig = analysis_utils.plot_monthly(\n", " _polygons_anomaly,\n", - " col='value',\n", - " cm = 'bwr',\n", - " vmin = -3, vmax = 3,\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS Anomaly [C]',\n", + " col=\"value\",\n", + " cm=\"bwr\",\n", + " vmin=-3,\n", + " vmax=3,\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS Anomaly [C]\",\n", ")\n", "\n", - "fig.savefig(\"2608_monthly_anomaly.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_monthly_anomaly.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -540,18 +590,20 @@ "metadata": {}, "outputs": [], "source": [ - "#6 month average anomaly\n", + "# 6 month average anomaly\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_anomaly,\n", - " col='value',\n", - " agg = 'mean',\n", - " cm = 'bwr',\n", - " vmin = -3.5, vmax = 3.5, n_colors = 7,\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS Anomaly [C]',\n", + " col=\"value\",\n", + " agg=\"mean\",\n", + " cm=\"bwr\",\n", + " vmin=-3.5,\n", + " vmax=3.5,\n", + " n_colors=7,\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS Anomaly [C]\",\n", ")\n", "\n", - "fig.savefig(\"2608_cum_anomaly.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_cum_anomaly.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -577,9 +629,15 @@ "metadata": {}, "outputs": [], "source": [ - "analysis_utils.make_csv(effect, socioeconomics, ensemble = True, \n", - " baseline_period=baseline_period, hotonly = \"coldonly\", \n", - " rate = False, age_weight = True)" + "analysis_utils.make_csv(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=\"coldonly\",\n", + " rate=False,\n", + " age_weight=True,\n", + ")" ] }, { @@ -599,60 +657,109 @@ "source": [ "def gdp_quartile(da_gdppc):\n", " # Determine quartiles from IR-level GDP per Capita (gdppc)\n", - " quartile_labels = pd.qcut(da_gdppc.values, q=4, labels=['low', 'mid-low', 'mid-high', 'high'])\n", - " quartile_da = xr.DataArray(\n", + " quartile_labels = pd.qcut(\n", + " da_gdppc.values, q=4, labels=[\"low\", \"mid-low\", \"mid-high\", \"high\"]\n", + " )\n", + " quartile_da = xr.DataArray(\n", " quartile_labels.astype(str),\n", " dims=da_gdppc.dims,\n", " coords=da_gdppc.coords,\n", - " name='gdp_quartile'\n", + " name=\"gdp_quartile\",\n", " )\n", - " return xr.Dataset({'gdppc': da_gdppc, 'gdp_quartile': quartile_da}) \n", - "\n", - "#Load age-weighted data (total deaths)\n", - "impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = True, baseline_period=baseline_period, hotonly = False, rate = False, age_weight = True)\n", - "_polygons_impact = analysis_utils.xarray_to_gpd(impact.mean(dim = 'number').sum(dim = 'month'), _polygons)\n", + " return xr.Dataset({\"gdppc\": da_gdppc, \"gdp_quartile\": quartile_da})\n", + "\n", + "\n", + "# Load age-weighted data (total deaths)\n", + "impact = analysis_utils.compute_impact(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=False,\n", + " rate=False,\n", + " age_weight=True,\n", + ")\n", + "_polygons_impact = analysis_utils.xarray_to_gpd(\n", + " impact.mean(dim=\"number\").sum(dim=\"month\"), _polygons\n", + ")\n", "# Population-weighted total deaths\n", - "df_country = _polygons_impact[['age_weighted_impact', 'ISO']].groupby('ISO').sum()\n", + "df_country = _polygons_impact[[\"age_weighted_impact\", \"ISO\"]].groupby(\"ISO\").sum()\n", + "\n", + "# Load data for all age cohorts\n", + "for cohort in [\"age0to4\", \"age5to64\", \"age65plus\"]:\n", + " impact = analysis_utils.compute_impact(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=False,\n", + " rate=False,\n", + " age_weight=False,\n", + " cohort=cohort,\n", + " )\n", + " _polygons_impact = analysis_utils.xarray_to_gpd(\n", + " impact.mean(dim=\"number\").sum(dim=\"month\"), _polygons\n", + " )\n", "\n", - "#Load data for all age cohorts\n", - "for cohort in ['age0to4', 'age5to64', 'age65plus']:\n", - " impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = True, baseline_period=baseline_period, hotonly = False, rate = False, age_weight = False, cohort = cohort)\n", - " _polygons_impact = analysis_utils.xarray_to_gpd(impact.mean(dim = 'number').sum(dim = 'month'), _polygons)\n", - " \n", - " col_i = f'{cohort}_impact'\n", - " df_country[col_i] = _polygons_impact[[col_i, 'ISO']].groupby('ISO').sum()\n", + " col_i = f\"{cohort}_impact\"\n", + " df_country[col_i] = _polygons_impact[[col_i, \"ISO\"]].groupby(\"ISO\").sum()\n", "\n", - "#World Bank Income Rankings\n", - "#income_class = pd.read_csv(os.path.join(DATA_DIR, \"income_class_2025.csv\")) \n", - "#out_df = df_country.merge(income_class, how = 'left', left_on = 'ISO', right_on = 'Code')\n", + "# World Bank Income Rankings\n", + "# income_class = pd.read_csv(os.path.join(DATA_DIR, \"income_class_2025.csv\"))\n", + "# out_df = df_country.merge(income_class, how = 'left', left_on = 'ISO', right_on = 'Code')\n", "\n", - "#Pop-weighted GDP by Country\n", + "# Pop-weighted GDP by Country\n", "socioeconomics = socioeconomics.load()\n", "gdp_quartile_df = gdp_quartile(\n", - " socioeconomics.groupby('iso3').map(lambda ds: ds['gdppc'].weighted(ds['pop']).mean())\n", + " socioeconomics.groupby(\"iso3\").map(\n", + " lambda ds: ds[\"gdppc\"].weighted(ds[\"pop\"]).mean()\n", + " )\n", ").to_dataframe()\n", "\n", - "pop_df = socioeconomics[['iso3', 'pop', 'pop0to4', 'pop5to64', 'pop65plus']].groupby('iso3').sum()\n", + "pop_df = (\n", + " socioeconomics[[\"iso3\", \"pop\", \"pop0to4\", \"pop5to64\", \"pop65plus\"]]\n", + " .groupby(\"iso3\")\n", + " .sum()\n", + ")\n", "\n", - "out_df = df_country.reset_index().merge(gdp_quartile_df.reset_index(), how = 'left', left_on = 'ISO', right_on = 'iso3')\n", - "out_df = out_df.merge(pop_df.to_dataframe().reset_index(), how = 'left', left_on = 'ISO', right_on = 'iso3')\n", + "out_df = df_country.reset_index().merge(\n", + " gdp_quartile_df.reset_index(), how=\"left\", left_on=\"ISO\", right_on=\"iso3\"\n", + ")\n", + "out_df = out_df.merge(\n", + " pop_df.to_dataframe().reset_index(), how=\"left\", left_on=\"ISO\", right_on=\"iso3\"\n", + ")\n", "\n", - "#Compute Mortality Rates from Total Mortality \n", - "out_df[\"age_weighted_impact_rate\"] = out_df['age_weighted_impact']*100000/out_df['pop']\n", - "for cohort in ['age0to4', 'age5to64', 'age65plus']:\n", - " col_t = f'{cohort}_impact_rate'\n", - " col_i = f'{cohort}_impact'\n", + "# Compute Mortality Rates from Total Mortality\n", + "out_df[\"age_weighted_impact_rate\"] = (\n", + " out_df[\"age_weighted_impact\"] * 100000 / out_df[\"pop\"]\n", + ")\n", + "for cohort in [\"age0to4\", \"age5to64\", \"age65plus\"]:\n", + " col_t = f\"{cohort}_impact_rate\"\n", + " col_i = f\"{cohort}_impact\"\n", " _cohortstem = cohort[3:]\n", " col_p = f\"pop{_cohortstem}\"\n", - " out_df[col_t] = out_df[col_i]*100000/out_df[col_p]\n", - "\n", - "out_df = out_df[[\"ISO\", \"gdppc\", \"gdp_quartile\",\n", - " \"age_weighted_impact\", \"age_weighted_impact_rate\", \"pop\",\n", - " \"age0to4_impact\", \"age0to4_impact_rate\", \"pop0to4\",\n", - " \"age5to64_impact\", \"age5to64_impact_rate\", \"pop5to64\",\n", - " \"age65plus_impact\", \"age65plus_impact_rate\", \"pop65plus\",\n", - " ]]\n", - "out_df.to_csv(\"2608_impact_country_gdppc.csv\", index = False)" + " out_df[col_t] = out_df[col_i] * 100000 / out_df[col_p]\n", + "\n", + "out_df = out_df[\n", + " [\n", + " \"ISO\",\n", + " \"gdppc\",\n", + " \"gdp_quartile\",\n", + " \"age_weighted_impact\",\n", + " \"age_weighted_impact_rate\",\n", + " \"pop\",\n", + " \"age0to4_impact\",\n", + " \"age0to4_impact_rate\",\n", + " \"pop0to4\",\n", + " \"age5to64_impact\",\n", + " \"age5to64_impact_rate\",\n", + " \"pop5to64\",\n", + " \"age65plus_impact\",\n", + " \"age65plus_impact_rate\",\n", + " \"pop65plus\",\n", + " ]\n", + "]\n", + "out_df.to_csv(\"2608_impact_country_gdppc.csv\", index=False)" ] }, { @@ -662,7 +769,7 @@ "metadata": {}, "outputs": [], "source": [ - "out_df.sort_values(\"age_weighted_impact_rate\", ascending = False)" + "out_df.sort_values(\"age_weighted_impact_rate\", ascending=False)" ] }, { @@ -682,10 +789,14 @@ "source": [ "### TEST DATA ###\n", "# IR Climate Data\n", - "forecast_gmfd = xr.open_zarr(os.path.join(DATA_DIR, \"daily_forecast_GMFD_tas.zarr\"))-273.15\n", - "reanalysis_gmfd = xr.open_zarr(os.path.join(DATA_DIR, \"daily_reanalysis_GMFD_tas.zarr\"))-273.15\n", + "forecast_gmfd = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"daily_forecast_GMFD_tas.zarr\")) - 273.15\n", + ")\n", + "reanalysis_gmfd = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"daily_reanalysis_GMFD_tas.zarr\")) - 273.15\n", + ")\n", "\n", - "#Intermediate Betas/MMT from Project.Py\n", + "# Intermediate Betas/MMT from Project.Py\n", "betas_mmt = xr.open_zarr(os.path.join(DATA_DIR, \"beta_mmt_GMFD.zarr\"))" ] }, @@ -696,21 +807,15 @@ "metadata": {}, "outputs": [], "source": [ - "import matplotlib.patheffects as pe\n", - "from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec\n", - "\n", - "import cartopy.crs as ccrs\n", - "from scipy.interpolate import interp1d\n", - "from shapely import box\n", - "\n", "import diagnostic_utils\n", + "\n", "importlib.reload(diagnostic_utils)\n", "\n", "### Analysis Options ###\n", - "rate = True # Mortality Rate vs Total Deaths\n", - "hotonly = False # Hot-Only Deaths\n", - "age_weight = False #Proportional Age-weighting vs single age group\n", - "density_hist = False #Controls daily data histograms (not computation normalization)\n", + "rate = True # Mortality Rate vs Total Deaths\n", + "hotonly = False # Hot-Only Deaths\n", + "age_weight = False # Proportional Age-weighting vs single age group\n", + "density_hist = False # Controls daily data histograms (not computation normalization)\n", "\n", "### Region Selection ###\n", "region_filters = [\n", @@ -732,27 +837,38 @@ "\n", "# Compute Impact\n", "projected = effect\n", - "regional_sum = analysis_utils.compute_impact(projected, socioeconomics, baseline_period= baseline_period, hotonly = hotonly, rate = rate, age_weight = age_weight)\n", + "regional_sum = analysis_utils.compute_impact(\n", + " projected,\n", + " socioeconomics,\n", + " baseline_period=baseline_period,\n", + " hotonly=hotonly,\n", + " rate=rate,\n", + " age_weight=age_weight,\n", + ")\n", "\n", - "#Get Months\n", + "# Get Months\n", "months = regional_sum.month.values\n", "\n", "### Load Climate Data ###\n", - "#Select representative ensemble member\n", - "forecast_slice = forecast_gmfd['tas'].sel(number = 1)\n", - "forecast_local = forecast_slice.sel(region = region_filter)\n", + "# Select representative ensemble member\n", + "forecast_slice = forecast_gmfd[\"tas\"].sel(number=1)\n", + "forecast_local = forecast_slice.sel(region=region_filter)\n", "\n", - "#Select representative year\n", - "reanalysis_slice = reanalysis_gmfd['tas'].sel(time=slice('2000-01-01', '2000-12-31'))\n", + "# Select representative year\n", + "reanalysis_slice = reanalysis_gmfd[\"tas\"].sel(time=slice(\"2000-01-01\", \"2000-12-31\"))\n", "month_mask = reanalysis_slice.month.isin(months).compute()\n", "reanalysis_slice = reanalysis_slice.where(month_mask, drop=True)\n", - "reanalysis_local = reanalysis_slice.sel(region = region_filter)\n", + "reanalysis_local = reanalysis_slice.sel(region=region_filter)\n", "\n", - "combined = np.concat([forecast_local.values.flatten(), reanalysis_local.values.flatten()])\n", + "combined = np.concat(\n", + " [forecast_local.values.flatten(), reanalysis_local.values.flatten()]\n", + ")\n", "xmin, xmax = np.percentile(combined, [1, 100])\n", "\n", "# Compute Cumulative Effect\n", - "cdf_data, max_cdf = diagnostic_utils.compute_cumulative_effect(forecast_slice, reanalysis_slice, region_filter, months, hotonly = hotonly)" + "cdf_data, max_cdf = diagnostic_utils.compute_cumulative_effect(\n", + " forecast_slice, reanalysis_slice, region_filter, months, hotonly=hotonly\n", + ")" ] }, { @@ -762,35 +878,63 @@ "metadata": {}, "outputs": [], "source": [ - "#Histogram\n", + "# Histogram\n", "\n", "b_color = \"#999999\"\n", "f_color = \"#ffbb6f\"\n", "r_color = \"#5e4c5f\"\n", "\n", "fig, (ax_hist, ax_cdf) = plt.subplots(\n", - " 2, 1, sharex=True, figsize=(6, 6),\n", - " gridspec_kw={'height_ratios': [3, 1], 'hspace': 0.08}\n", + " 2,\n", + " 1,\n", + " sharex=True,\n", + " figsize=(6, 6),\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"hspace\": 0.08},\n", ")\n", "\n", "\n", - "ref_vals = betas_mmt['mmt'].sel(region=region_filter).sel(age_cohort='age65plus')\n", - "betas = betas_mmt['beta_hotonly'] if hotonly else betas_mmt['beta']\n", - "da_temp_bins = betas.sel(region=region_filter).sel(age_cohort='age65plus')\n", - "da_temp_bins['tas_bin'] = da_temp_bins['tas_bin'] \n", + "ref_vals = betas_mmt[\"mmt\"].sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + "betas = betas_mmt[\"beta_hotonly\"] if hotonly else betas_mmt[\"beta\"]\n", + "da_temp_bins = betas.sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + "da_temp_bins[\"tas_bin\"] = da_temp_bins[\"tas_bin\"]\n", "\n", "# Density histograms on hidden twin axis\n", - "bins = da_temp_bins.coords['tas_bin']\n", - "ax_hist.hist(reanalysis_local.values.flatten(), bins=bins, alpha=0.5, label=\"ERA5 (2000)\", density=density_hist, color=b_color)\n", - "ax_hist.hist(forecast_local.values.flatten(), bins=bins, alpha=0.5, label=\"Forecast\", density=density_hist, color=f_color)\n", - "ax_hist.axvline(forecast_local.values.flatten().mean(), color=f_color, linewidth=1, linestyle='--')\n", - "ax_hist.axvline(reanalysis_local.values.flatten().mean(), color=b_color, linewidth=1, linestyle='--')\n", + "bins = da_temp_bins.coords[\"tas_bin\"]\n", + "ax_hist.hist(\n", + " reanalysis_local.values.flatten(),\n", + " bins=bins,\n", + " alpha=0.5,\n", + " label=\"ERA5 (2000)\",\n", + " density=density_hist,\n", + " color=b_color,\n", + ")\n", + "ax_hist.hist(\n", + " forecast_local.values.flatten(),\n", + " bins=bins,\n", + " alpha=0.5,\n", + " label=\"Forecast\",\n", + " density=density_hist,\n", + " color=f_color,\n", + ")\n", + "ax_hist.axvline(\n", + " forecast_local.values.flatten().mean(), color=f_color, linewidth=1, linestyle=\"--\"\n", + ")\n", + "ax_hist.axvline(\n", + " reanalysis_local.values.flatten().mean(), color=b_color, linewidth=1, linestyle=\"--\"\n", + ")\n", "ax_hist.set_ylabel(\"Count [Days]\")\n", "ax_hist.set_xlim(xmin, xmax)\n", "\n", "# Response curve on the primary (left) axis\n", "ax_hist2 = ax_hist.twinx()\n", - "ax_hist2.step(da_temp_bins.coords['tas_bin'], da_temp_bins.values, where='pre', color=r_color, linewidth=2, label=\"Response curve\")\n", + "ax_hist2.step(\n", + " da_temp_bins.coords[\"tas_bin\"],\n", + " da_temp_bins.values,\n", + " where=\"pre\",\n", + " color=r_color,\n", + " linewidth=2,\n", + " label=\"Response curve\",\n", + ")\n", "ax_hist2.axvline(ref_vals.values, color=r_color, linewidth=1, label=\"MMT\")\n", "ax_hist2.set_ylabel(\"Daily Response 65+ \\n [deaths/100k]\")\n", "ax_hist2.set_ylim(0, 10)\n", @@ -800,23 +944,36 @@ "era5_centers, era5_cdf, fc_centers, fc_cdf = cdf_data\n", "\n", "ax_cdf.set_ylabel(\"6-Month Cum.\\n [Deaths/100k]\")\n", - "ax_cdf.set_ylim(0, max_cdf+2)\n", + "ax_cdf.set_ylim(0, max_cdf + 2)\n", "ax_cdf.set_xlim(xmin, xmax)\n", "ax_cdf.set_xlabel(\"Tas [C]\")\n", - "ax_cdf.step(fc_centers, fc_cdf, where='post', color=f_color, label=\"Forecast CDF\")\n", - "ax_cdf.step(era5_centers, era5_cdf, where='post', color=b_color, label=\"ERA5 CDF\")\n", - "\n", - "#Legend\n", - "handles = (ax_hist.get_legend_handles_labels()[0]\n", - " + ax_hist2.get_legend_handles_labels()[0]\n", - " + ax_cdf.get_legend_handles_labels()[0])\n", - "labels = (ax_hist.get_legend_handles_labels()[1]\n", - " + ax_hist2.get_legend_handles_labels()[1]\n", - " + ax_cdf.get_legend_handles_labels()[1])\n", - "ax_cdf.legend(handles, labels, fontsize=7, loc='upper center',\n", - " bbox_to_anchor=(0.5, -0.35), ncol=3, frameon=False)\n", - "ax_hist.set_title(f'{region_filter}')\n", - "fig.savefig(f'response_diagnostic_{region_filter}_2607.png', dpi=600, bbox_inches=\"tight\")" + "ax_cdf.step(fc_centers, fc_cdf, where=\"post\", color=f_color, label=\"Forecast CDF\")\n", + "ax_cdf.step(era5_centers, era5_cdf, where=\"post\", color=b_color, label=\"ERA5 CDF\")\n", + "\n", + "# Legend\n", + "handles = (\n", + " ax_hist.get_legend_handles_labels()[0]\n", + " + ax_hist2.get_legend_handles_labels()[0]\n", + " + ax_cdf.get_legend_handles_labels()[0]\n", + ")\n", + "labels = (\n", + " ax_hist.get_legend_handles_labels()[1]\n", + " + ax_hist2.get_legend_handles_labels()[1]\n", + " + ax_cdf.get_legend_handles_labels()[1]\n", + ")\n", + "ax_cdf.legend(\n", + " handles,\n", + " labels,\n", + " fontsize=7,\n", + " loc=\"upper center\",\n", + " bbox_to_anchor=(0.5, -0.35),\n", + " ncol=3,\n", + " frameon=False,\n", + ")\n", + "ax_hist.set_title(f\"{region_filter}\")\n", + "fig.savefig(\n", + " f\"response_diagnostic_{region_filter}_2607.png\", dpi=600, bbox_inches=\"tight\"\n", + ")" ] }, { diff --git a/analysis/a03_systems_comp.ipynb b/analysis/a03_systems_comp.ipynb index da21a58..4c4f34b 100644 --- a/analysis/a03_systems_comp.ipynb +++ b/analysis/a03_systems_comp.ipynb @@ -19,17 +19,19 @@ "metadata": {}, "outputs": [], "source": [ - "centres = {'ecmwf': ['51'], \n", - " \"meteo_france\": ['9'], \n", - " 'ukmo': ['610'], \n", - " 'dwd': ['22'], \n", - " 'cmcc': ['4'], \n", - " 'ncep': ['2'], \n", - " 'jma': ['4'], \n", - " 'eccc': ['4','5'],\n", - " 'bom': ['2']}\n", + "centres = {\n", + " \"ecmwf\": [\"51\"],\n", + " \"meteo_france\": [\"9\"],\n", + " \"ukmo\": [\"610\"],\n", + " \"dwd\": [\"22\"],\n", + " \"cmcc\": [\"4\"],\n", + " \"ncep\": [\"2\"],\n", + " \"jma\": [\"4\"],\n", + " \"eccc\": [\"4\", \"5\"],\n", + " \"bom\": [\"2\"],\n", + "}\n", "\n", - "region_sel = 'SAS'" + "region_sel = \"SAS\"" ] }, { @@ -50,6 +52,7 @@ " \"fcstdev\": \"monthly_standard_deviation\",\n", "}\n", "\n", + "\n", "def zip_to_datatree(zip_path, extract_dir):\n", " with zipfile.ZipFile(zip_path) as z:\n", " names = z.namelist()\n", @@ -64,6 +67,7 @@ "\n", " return xr.DataTree.from_dict(tree_dict)\n", "\n", + "\n", "def land_only(da):\n", " land_mask = regionmask.defined_regions.natural_earth_v5_0_0.land_110.mask(da)\n", " da_land_only = da.where(land_mask.notnull())\n", @@ -80,10 +84,10 @@ "trees = {}\n", "for centre, systems in centres.items():\n", " for system in systems:\n", - " path = f'/home/emily_zuetell/projects/poreallas/scripts//data/raw/2607_{centre}_s{system}.nc'\n", + " path = f\"/home/emily_zuetell/projects/poreallas/scripts//data/raw/2607_{centre}_s{system}.nc\"\n", " extract_dir = path.replace(\".nc\", \"_extracted\")\n", "\n", - " var = f'{centre}_s{system}'\n", + " var = f\"{centre}_s{system}\"\n", " trees[var] = zip_to_datatree(path, extract_dir)" ] }, @@ -102,9 +106,11 @@ "metadata": {}, "outputs": [], "source": [ - "#TODO switch to xesmf\n", + "# TODO switch to xesmf\n", "def regrid_to_target(da, target):\n", - " return da.interp(latitude=target.latitude, longitude=target.longitude, method=\"linear\")\n" + " return da.interp(\n", + " latitude=target.latitude, longitude=target.longitude, method=\"linear\"\n", + " )" ] }, { @@ -148,7 +154,7 @@ "metadata": {}, "outputs": [], "source": [ - "diff = combined_anomaly.sel(centre_system = 'ecmwf_s51')-combined_anomaly" + "diff = combined_anomaly.sel(centre_system=\"ecmwf_s51\") - combined_anomaly" ] }, { @@ -170,9 +176,8 @@ ], "source": [ "import cartopy.crs as ccrs\n", - "import matplotlib.pyplot as plt\n", "\n", - "p = diff.sel(forecastMonth = 2).plot(\n", + "p = diff.sel(forecastMonth=2).plot(\n", " col=\"centre_system\",\n", " col_wrap=3,\n", " vmin=-3,\n", @@ -182,7 +187,7 @@ " subplot_kws={\"projection\": ccrs.PlateCarree()},\n", " transform=ccrs.PlateCarree(),\n", ")\n", - "p.fig.suptitle('ForecastMonth = 2')\n", + "p.fig.suptitle(\"ForecastMonth = 2\")\n", "for ax in p.axs.flat:\n", " ax.coastlines()" ] @@ -202,7 +207,7 @@ "metadata": {}, "outputs": [], "source": [ - "system_std = combined_anomaly.std('centre_system')" + "system_std = combined_anomaly.std(\"centre_system\")" ] }, { @@ -293,24 +298,25 @@ "outputs": [], "source": [ "import pandas as pd\n", - "import regionmask\n", "\n", "# e.g. AR6 land regions\n", "ar6_regions = regionmask.defined_regions.ar6.land\n", "\n", + "\n", "def select_region(da):\n", " mask = ar6_regions.mask(da)\n", " region_da = da.where(mask == ar6_regions.map_keys(region_sel))\n", " return region_da\n", "\n", + "\n", "records = []\n", "for centre_system, tree in trees.items(): # e.g. key = \"ecmwf_s51\"\n", - " da = tree[\"monthly_mean\"][\"t2m\"] \n", + " da = tree[\"monthly_mean\"][\"t2m\"]\n", " da_hc = tree[\"hindcast_climate_mean\"][\"t2m\"]\n", - " #select region\n", + " # select region\n", " da = select_region(da)\n", " da_hc = select_region(da_hc)\n", - " df = (da-da_hc).to_dataframe().reset_index()\n", + " df = (da - da_hc).to_dataframe().reset_index()\n", " df[\"centre_system\"] = centre_system\n", " records.append(df[[\"centre_system\", \"forecastMonth\", \"t2m\"]])\n", "\n", @@ -328,11 +334,11 @@ "for centre_system, tree in trees.items():\n", " da = tree[\"ensemble_mean\"][\"t2m\"]\n", " da_hc = tree[\"hindcast_climate_mean\"][\"t2m\"]\n", - " #select region\n", + " # select region\n", " da = select_region(da)\n", " da_hc = select_region(da_hc)\n", "\n", - " df = (da-da_hc).to_dataframe().reset_index()\n", + " df = (da - da_hc).to_dataframe().reset_index()\n", " df[\"centre_system\"] = centre_system\n", " records_em.append(df[[\"centre_system\", \"forecastMonth\", \"t2m\"]])\n", "\n", @@ -360,11 +366,11 @@ ")\n", "for ax, month in zip(g.axes.flat, g.col_names):\n", " vals = em_df.loc[em_df[\"forecastMonth\"] == month, \"t2m\"]\n", - " #q1, med, q3 = vals.quantile([0.25, 0.5, 0.75])\n", - " #ax.axvspan(q1, q3, color=\"black\", alpha=0.15, zorder=10)\n", + " # q1, med, q3 = vals.quantile([0.25, 0.5, 0.75])\n", + " # ax.axvspan(q1, q3, color=\"black\", alpha=0.15, zorder=10)\n", " ax.axvline(vals.mean(), color=\"black\", linewidth=1.5, zorder=11)\n", - " ax.set_xlabel('t2m anomaly')\n", - "g.fig.suptitle(f'T2M Anomaly: AR6 Region = {region_sel}', y=1.02)\n", + " ax.set_xlabel(\"t2m anomaly\")\n", + "g.fig.suptitle(f\"T2M Anomaly: AR6 Region = {region_sel}\", y=1.02)\n", "g.fig" ] }, diff --git a/analysis/a06_postcard_map.ipynb b/analysis/a06_postcard_map.ipynb index 18e840c..f89d6c9 100644 --- a/analysis/a06_postcard_map.ipynb +++ b/analysis/a06_postcard_map.ipynb @@ -45,11 +45,11 @@ " \"experiment\": \"ssp2_4_5\",\n", " \"domain\": \"global\",\n", " \"period\": \"2015-2100\",\n", - " \"variable\": \"monthly_mean_of_daily_mean_temperature\"\n", + " \"variable\": \"monthly_mean_of_daily_mean_temperature\",\n", "}\n", "\n", "client = cdsapi.Client()\n", - "client.retrieve(dataset, request, \"./data/cmip6_2015_2100_monthly_tmean_ssp245.nc\")\n" + "client.retrieve(dataset, request, \"./data/cmip6_2015_2100_monthly_tmean_ssp245.nc\")" ] }, { @@ -66,11 +66,11 @@ " \"experiment\": \"historical\",\n", " \"domain\": \"global\",\n", " \"period\": \"1850-2014\",\n", - " \"variable\": \"monthly_mean_of_daily_mean_temperature\"\n", + " \"variable\": \"monthly_mean_of_daily_mean_temperature\",\n", "}\n", "\n", "client = cdsapi.Client()\n", - "client.retrieve(dataset, request, \"./data/cmip6_hist_1850_2014_monthly_tmean.nc\")\n" + "client.retrieve(dataset, request, \"./data/cmip6_hist_1850_2014_monthly_tmean.nc\")" ] }, { @@ -88,37 +88,50 @@ " \"originating_centre\": \"ecmwf\",\n", " \"system\": \"51\",\n", " \"variable\": [\"2m_temperature\"],\n", - " \"product_type\": [\n", - " \"monthly_mean\"\n", - " ],\n", + " \"product_type\": [\"monthly_mean\"],\n", " \"year\": [\n", - " \"1981\", \"1982\", \"1983\",\n", - " \"1984\", \"1985\", \"1986\",\n", - " \"1987\", \"1988\", \"1989\",\n", - " \"1990\", \"1991\", \"1992\",\n", - " \"1993\", \"1994\", \"1995\",\n", - " \"1996\", \"1997\", \"1998\",\n", - " \"1999\", \"2000\", \"2001\",\n", - " \"2002\", \"2003\", \"2004\",\n", - " \"2005\", \"2006\", \"2007\",\n", - " \"2008\", \"2009\", \"2010\",\n", - " \"2011\", \"2012\", \"2013\",\n", - " \"2014\"\n", + " \"1981\",\n", + " \"1982\",\n", + " \"1983\",\n", + " \"1984\",\n", + " \"1985\",\n", + " \"1986\",\n", + " \"1987\",\n", + " \"1988\",\n", + " \"1989\",\n", + " \"1990\",\n", + " \"1991\",\n", + " \"1992\",\n", + " \"1993\",\n", + " \"1994\",\n", + " \"1995\",\n", + " \"1996\",\n", + " \"1997\",\n", + " \"1998\",\n", + " \"1999\",\n", + " \"2000\",\n", + " \"2001\",\n", + " \"2002\",\n", + " \"2003\",\n", + " \"2004\",\n", + " \"2005\",\n", + " \"2006\",\n", + " \"2007\",\n", + " \"2008\",\n", + " \"2009\",\n", + " \"2010\",\n", + " \"2011\",\n", + " \"2012\",\n", + " \"2013\",\n", + " \"2014\",\n", " ],\n", " \"month\": [\"07\"],\n", - " \"leadtime_month\": [\n", - " \"1\",\n", - " \"2\",\n", - " \"3\",\n", - " \"4\",\n", - " \"5\",\n", - " \"6\"\n", - " ],\n", - " \"data_format\": \"netcdf\"\n", + " \"leadtime_month\": [\"1\", \"2\", \"3\", \"4\", \"5\", \"6\"],\n", + " \"data_format\": \"netcdf\",\n", "}\n", "\n", "client = cdsapi.Client()\n", - "client.retrieve(dataset, request, \"./data/seas5_1981_2014_monthly_tmean_2.nc\")\n" + "client.retrieve(dataset, request, \"./data/seas5_1981_2014_monthly_tmean_2.nc\")" ] }, { @@ -134,22 +147,11 @@ " \"originating_centre\": \"ecmwf\",\n", " \"system\": \"51\",\n", " \"variable\": [\"2m_temperature\"],\n", - " \"product_type\": [\n", - " \"monthly_mean\"\n", - " ],\n", - " \"year\": [\n", - " \"2026\"\n", - " ],\n", + " \"product_type\": [\"monthly_mean\"],\n", + " \"year\": [\"2026\"],\n", " \"month\": [\"07\"],\n", - " \"leadtime_month\": [\n", - " \"1\",\n", - " \"2\",\n", - " \"3\",\n", - " \"4\",\n", - " \"5\",\n", - " \"6\"\n", - " ],\n", - " \"data_format\": \"netcdf\"\n", + " \"leadtime_month\": [\"1\", \"2\", \"3\", \"4\", \"5\", \"6\"],\n", + " \"data_format\": \"netcdf\",\n", "}\n", "\n", "client = cdsapi.Client()\n", @@ -172,6 +174,7 @@ "outputs": [], "source": [ "import zipfile\n", + "\n", "with zipfile.ZipFile(\"./data/seas5_1981_2014_monthly_tmean_2.nc\") as z:\n", " z.extractall(\"./data/\")" ] @@ -202,7 +205,7 @@ "metadata": {}, "outputs": [], "source": [ - "#Define Forecast Months\n", + "# Define Forecast Months\n", "fc_months = [8, 9, 10, 11, 12, 1]" ] }, @@ -233,8 +236,14 @@ ], "source": [ "## CMIP6 Projections\n", - "ar6_project = xr.open_dataset(os.path.join(DATA_DIR,\"t_CMIP6_ssp370_mon_201501-210012.nc\"), chunks={\"member\": -1, \"time\": -1, \"lat\": 30, \"lon\": 30})\n", - "ar6_hist = xr.open_dataset(os.path.join(DATA_DIR,\"t_CMIP6_historical_mon_185001-201412.nc\"), chunks={\"member\": -1, \"time\": -1, \"lat\": 30, \"lon\": 30})" + "ar6_project = xr.open_dataset(\n", + " os.path.join(DATA_DIR, \"t_CMIP6_ssp370_mon_201501-210012.nc\"),\n", + " chunks={\"member\": -1, \"time\": -1, \"lat\": 30, \"lon\": 30},\n", + ")\n", + "ar6_hist = xr.open_dataset(\n", + " os.path.join(DATA_DIR, \"t_CMIP6_historical_mon_185001-201412.nc\"),\n", + " chunks={\"member\": -1, \"time\": -1, \"lat\": 30, \"lon\": 30},\n", + ")" ] }, { @@ -244,9 +253,14 @@ "metadata": {}, "outputs": [], "source": [ - "#IPCC Recent Baseline Period 1995-2014\n", - "hist_mean = ar6_hist['t'].sel(time = slice(\"1995-01-01\", \"2014-12-31\")).groupby(\"time.month\").mean()\n", - "hist_mean = hist_mean.sel(month = fc_months).mean(dim = 'member').compute()" + "# IPCC Recent Baseline Period 1995-2014\n", + "hist_mean = (\n", + " ar6_hist[\"t\"]\n", + " .sel(time=slice(\"1995-01-01\", \"2014-12-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + ")\n", + "hist_mean = hist_mean.sel(month=fc_months).mean(dim=\"member\").compute()" ] }, { @@ -256,10 +270,12 @@ "metadata": {}, "outputs": [], "source": [ - "projection = ar6_project['t']\n", - "projection = projection.assign_coords(year=projection.time.dt.year, month=projection.time.dt.month)\n", + "projection = ar6_project[\"t\"]\n", + "projection = projection.assign_coords(\n", + " year=projection.time.dt.year, month=projection.time.dt.month\n", + ")\n", "projection = projection.set_index(time=[\"year\", \"month\"]).unstack(\"time\")\n", - "projection = projection.sel(month = fc_months).mean(dim = ['member']).compute()" + "projection = projection.sel(month=fc_months).mean(dim=[\"member\"]).compute()" ] }, { @@ -270,10 +286,24 @@ "outputs": [], "source": [ "### Seasonal Forecast ###\n", - "#Raw Forecast/Hindcast\n", - "s51 = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/parsed/08_ecmwf_parsed.zarr\")\n", - "forecast = s51.sel(time=slice('2026-08-01', '2027-01-31')).groupby(\"time.month\").mean().compute()\n", - "hindcast = s51.sel(time=slice('1995-08-01', '2014-01-31')).groupby(\"time.month\").mean().sel(month = [1, 8, 9, 10, 11, 12]).mean(dim = 'number').compute()\n" + "# Raw Forecast/Hindcast\n", + "s51 = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/parsed/08_ecmwf_parsed.zarr\"\n", + ")\n", + "forecast = (\n", + " s51.sel(time=slice(\"2026-08-01\", \"2027-01-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .compute()\n", + ")\n", + "hindcast = (\n", + " s51.sel(time=slice(\"1995-08-01\", \"2014-01-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[1, 8, 9, 10, 11, 12])\n", + " .mean(dim=\"number\")\n", + " .compute()\n", + ")" ] }, { @@ -293,30 +323,30 @@ "source": [ "## Statistics ##\n", "# Climate Projection\n", - "mean_proj = projection.mean(dim = 'member')\n", - "var_proj = (projection.std(dim = 'member'))**2\n", + "mean_proj = projection.mean(dim=\"member\")\n", + "var_proj = (projection.std(dim=\"member\")) ** 2\n", "n_proj = len(projection.member)\n", "\n", - "mean_base = hist_mean.mean(dim = 'member')\n", - "var_base = (hist_mean.std(dim = 'member'))**2\n", + "mean_base = hist_mean.mean(dim=\"member\")\n", + "var_base = (hist_mean.std(dim=\"member\")) ** 2\n", "n_base = len(hist_mean.member)\n", "\n", - "anomaly_proj = mean_proj-mean_base\n", - "#Rolling avg\n", + "anomaly_proj = mean_proj - mean_base\n", + "# Rolling avg\n", "anomaly_proj = anomaly_proj.rolling(year=5, center=False, min_periods=5).mean()\n", - "var_anomaly_proj = (var_proj/n_proj)+(var_base/n_base)\n", + "var_anomaly_proj = (var_proj / n_proj) + (var_base / n_base)\n", "\n", "# Seasonal Forecast\n", - "mean_fc = forecast['t2m'].mean(dim = 'number')\n", - "var_fc = (forecast['t2m'].var(dim = 'number'))\n", - "n_fc = len(forecast['t2m'].number)\n", + "mean_fc = forecast[\"t2m\"].mean(dim=\"number\")\n", + "var_fc = forecast[\"t2m\"].var(dim=\"number\")\n", + "n_fc = len(forecast[\"t2m\"].number)\n", "\n", - "mean_hc = hindcast['t2m'].mean(dim = 'number')\n", - "var_hc = (hindcast['t2m'].var(dim = 'number'))\n", - "n_hc = len(hindcast['t2m'].number)\n", + "mean_hc = hindcast[\"t2m\"].mean(dim=\"number\")\n", + "var_hc = hindcast[\"t2m\"].var(dim=\"number\")\n", + "n_hc = len(hindcast[\"t2m\"].number)\n", "\n", - "anomaly_fc = mean_fc-mean_hc\n", - "var_anomaly_fc = (var_fc/n_fc)+(var_hc/n_hc)" + "anomaly_fc = mean_fc - mean_hc\n", + "var_anomaly_fc = (var_fc / n_fc) + (var_hc / n_hc)" ] }, { @@ -327,10 +357,10 @@ "outputs": [], "source": [ "## Probability of Surpassing Projection Anomaly\n", - "diff = anomaly_fc-anomaly_proj\n", + "diff = anomaly_fc - anomaly_proj\n", "var_diff = var_anomaly_fc + var_anomaly_proj\n", - "#var_diff = var_anomaly_fc #Forecast uncertainty only (No projection uncertainty)\n", - "z = diff/np.sqrt(var_diff)\n", + "# var_diff = var_anomaly_fc #Forecast uncertainty only (No projection uncertainty)\n", + "z = diff / np.sqrt(var_diff)\n", "p_surpass = xr.apply_ufunc(norm.cdf, z)" ] }, @@ -362,15 +392,33 @@ } ], "source": [ - "#lat, lon = 42, 49\n", + "# lat, lon = 42, 49\n", "lat, lon = 9.9312, 76.2673\n", "fig, ax = plt.subplots()\n", - "p_surpass.sel(lat=lat, lon = lon, method = 'nearest').plot(ax = ax, label = \"p(fc>proj)\")\n", - "diff.sel(lat=lat, lon = lon, method = 'nearest').plot(ax = ax, label = \"diff(anomaly)\")\n", - "ax.axhline(0, color = 'k')\n", - "ax.axvline(p_surpass.sel(lat=lat, lon = lon, method = 'nearest').where(p_surpass.sel(lat=lat, lon = lon, method = 'nearest') >= 0.5).idxmin('year'), color = 'red', label = 'crossing year')\n", - "ax.axvline(p_surpass.sel(lat=lat, lon = lon, method = 'nearest').where(p_surpass.sel(lat=lat, lon = lon, method = 'nearest') >= 0.025).idxmin('year'), color = 'red', linestyle = '--')\n", - "ax.axvline(p_surpass.sel(lat=lat, lon = lon, method = 'nearest').where(p_surpass.sel(lat=lat, lon = lon, method = 'nearest') >= 0.975).idxmin('year'), color = 'red', linestyle = '--')\n", + "p_surpass.sel(lat=lat, lon=lon, method=\"nearest\").plot(ax=ax, label=\"p(fc>proj)\")\n", + "diff.sel(lat=lat, lon=lon, method=\"nearest\").plot(ax=ax, label=\"diff(anomaly)\")\n", + "ax.axhline(0, color=\"k\")\n", + "ax.axvline(\n", + " p_surpass.sel(lat=lat, lon=lon, method=\"nearest\")\n", + " .where(p_surpass.sel(lat=lat, lon=lon, method=\"nearest\") >= 0.5)\n", + " .idxmin(\"year\"),\n", + " color=\"red\",\n", + " label=\"crossing year\",\n", + ")\n", + "ax.axvline(\n", + " p_surpass.sel(lat=lat, lon=lon, method=\"nearest\")\n", + " .where(p_surpass.sel(lat=lat, lon=lon, method=\"nearest\") >= 0.025)\n", + " .idxmin(\"year\"),\n", + " color=\"red\",\n", + " linestyle=\"--\",\n", + ")\n", + "ax.axvline(\n", + " p_surpass.sel(lat=lat, lon=lon, method=\"nearest\")\n", + " .where(p_surpass.sel(lat=lat, lon=lon, method=\"nearest\") >= 0.975)\n", + " .idxmin(\"year\"),\n", + " color=\"red\",\n", + " linestyle=\"--\",\n", + ")\n", "ax.legend()" ] }, @@ -391,21 +439,23 @@ "source": [ "## Statistics ##\n", "# Climate Projection\n", - "mean_proj = projection.mean(dim = 'month')\n", - "mean_base = hist_mean.mean(dim = 'month')\n", + "mean_proj = projection.mean(dim=\"month\")\n", + "mean_base = hist_mean.mean(dim=\"month\")\n", "\n", - "anomaly_proj = mean_proj-mean_base\n", - "#Rolling avg\n", - "#anomaly_proj = anomaly_proj.rolling(year=8, center=False, min_periods=5).mean()\n", - "#Enforce monotonically increasing projection\n", + "anomaly_proj = mean_proj - mean_base\n", + "# Rolling avg\n", + "# anomaly_proj = anomaly_proj.rolling(year=8, center=False, min_periods=5).mean()\n", + "# Enforce monotonically increasing projection\n", "anomaly_proj_smooth = xr.apply_ufunc(\n", - " np.maximum.accumulate, anomaly_proj,\n", - " input_core_dims=[['year']], output_core_dims=[['year']],\n", - " kwargs={'axis': -1},\n", + " np.maximum.accumulate,\n", + " anomaly_proj,\n", + " input_core_dims=[[\"year\"]],\n", + " output_core_dims=[[\"year\"]],\n", + " kwargs={\"axis\": -1},\n", ")\n", "\n", "ds1 = anomaly_proj_smooth\n", - "decade1 = ((ds1['year'] // 10) * 10).rename('decade')\n", + "decade1 = ((ds1[\"year\"] // 10) * 10).rename(\"decade\")\n", "decade_min = ds1.groupby(decade1).min()\n", "decade_max = ds1.groupby(decade1).max()" ] @@ -438,10 +488,10 @@ } ], "source": [ - "anomaly_proj_smooth.isel(lat = 5, lon = 215).plot()\n", - "plt.axhline(ds2.isel(lat = 5, lon = 215).values)\n", - "plt.axvline(2040, color = 'black')\n", - "plt.axvline(2050, color = 'black')" + "anomaly_proj_smooth.isel(lat=5, lon=215).plot()\n", + "plt.axhline(ds2.isel(lat=5, lon=215).values)\n", + "plt.axvline(2040, color=\"black\")\n", + "plt.axvline(2050, color=\"black\")" ] }, { @@ -451,17 +501,17 @@ "metadata": {}, "outputs": [], "source": [ - "ds2 = forecast['tas']-hindcast['tas']\n", - "ds2 = ds2.mean(dim = 'month')\n", - "ds2 = isku_utils.lon_adjust(ds2, roll = True)\n", + "ds2 = forecast[\"tas\"] - hindcast[\"tas\"]\n", + "ds2 = ds2.mean(dim=\"month\")\n", + "ds2 = isku_utils.lon_adjust(ds2, roll=True)\n", "\n", - "#Add a delta \n", + "# Add a delta\n", "delta = 0.2\n", - "in_range = (ds2 >= (decade_min-delta)) & (ds2 <= (decade_max+delta))\n", - "below = ds2 < decade_min.min('decade')\n", - "above = ds2 > decade_max.max('decade')\n", + "in_range = (ds2 >= (decade_min - delta)) & (ds2 <= (decade_max + delta))\n", + "below = ds2 < decade_min.min(\"decade\")\n", + "above = ds2 > decade_max.max(\"decade\")\n", "\n", - "mode_decade = in_range.where(in_range ==True).idxmin('decade').compute()\n", + "mode_decade = in_range.where(in_range == True).idxmin(\"decade\").compute()\n", "\n", "da = xr.where(below, 2010, xr.where(above, 2110, mode_decade)).compute()" ] @@ -475,19 +525,19 @@ "source": [ "from scipy import stats\n", "\n", - "#Ensemble Stats\n", - "median = da.median(dim = 'number')\n", + "# Ensemble Stats\n", + "median = da.median(dim=\"number\")\n", "\n", "mode = xr.apply_ufunc(\n", " lambda x: stats.mode(x, axis=-1, keepdims=False).mode,\n", " da,\n", - " input_core_dims=[['number']],\n", + " input_core_dims=[[\"number\"]],\n", " vectorize=True,\n", ")\n", "mode_count = xr.apply_ufunc(\n", " lambda x: stats.mode(x, axis=-1, keepdims=False).count,\n", " da,\n", - " input_core_dims=[['number']],\n", + " input_core_dims=[[\"number\"]],\n", " vectorize=True,\n", ")" ] @@ -507,34 +557,36 @@ "metadata": {}, "outputs": [], "source": [ - "#Empirical Version\n", + "# Empirical Version\n", "# Climate Projection\n", - "mean_proj = projection.mean(dim = 'member')\n", - "mean_base = hist_mean.mean(dim = 'member')\n", - "anomaly_proj = mean_proj-mean_base\n", + "mean_proj = projection.mean(dim=\"member\")\n", + "mean_base = hist_mean.mean(dim=\"member\")\n", + "anomaly_proj = mean_proj - mean_base\n", "\n", "anomaly_proj = anomaly_proj.rolling(year=5, center=False, min_periods=5).mean()\n", - "emp_diff = (forecast['t2m']-hindcast['t2m']-anomaly_proj)\n", - "years_emp = emp_diff.where(emp_diff <= 0).idxmax('year')\n", + "emp_diff = forecast[\"t2m\"] - hindcast[\"t2m\"] - anomaly_proj\n", + "years_emp = emp_diff.where(emp_diff <= 0).idxmax(\"year\")\n", "\n", - "#Save Out\n", - "quantiles = years_emp.quantile([0.1, 0.5, 0.9], dim='number')\n", - "crossing_year = quantiles.sel(quantile = 0.5).drop_vars('quantile')\n", - "ci_low = quantiles.sel(quantile = 0.1).drop_vars('quantile')\n", - "ci_high = quantiles.sel(quantile = 0.9).drop_vars('quantile')\n", + "# Save Out\n", + "quantiles = years_emp.quantile([0.1, 0.5, 0.9], dim=\"number\")\n", + "crossing_year = quantiles.sel(quantile=0.5).drop_vars(\"quantile\")\n", + "ci_low = quantiles.sel(quantile=0.1).drop_vars(\"quantile\")\n", + "ci_high = quantiles.sel(quantile=0.9).drop_vars(\"quantile\")\n", "ci_low = ci_low.fillna(2014)\n", "ci_high = ci_high.fillna(2101)\n", "\n", "ci_width = ci_high - ci_low\n", "\n", - "ds_out = xr.Dataset({\n", - " 'crossing_year': crossing_year,\n", - " 'ci_width_95': ci_width,\n", - " 'ci_10': ci_low,\n", - " 'ci_90':ci_high\n", - "})\n", + "ds_out = xr.Dataset(\n", + " {\n", + " \"crossing_year\": crossing_year,\n", + " \"ci_width_95\": ci_width,\n", + " \"ci_10\": ci_low,\n", + " \"ci_90\": ci_high,\n", + " }\n", + ")\n", "\n", - "#ds_out.to_zarr('./data/2607_fc_anomaly_year_fcvarOnly_emp_ssp245.zarr')" + "# ds_out.to_zarr('./data/2607_fc_anomaly_year_fcvarOnly_emp_ssp245.zarr')" ] }, { @@ -565,24 +617,36 @@ } ], "source": [ - "#lat, lon = 42, 49\n", + "# lat, lon = 42, 49\n", "lat, lon = 9.9312, 76.2673\n", "fig, ax = plt.subplots()\n", - "data = (forecast['t2m']-hindcast['t2m']).sel(lat=lat, lon = lon, method = 'nearest')\n", - "anomaly_proj.sel(lat=lat, lon = lon, method = 'nearest').plot(ax = ax, label = \"anomaly_proj\")\n", + "data = (forecast[\"t2m\"] - hindcast[\"t2m\"]).sel(lat=lat, lon=lon, method=\"nearest\")\n", + "anomaly_proj.sel(lat=lat, lon=lon, method=\"nearest\").plot(ax=ax, label=\"anomaly_proj\")\n", + "\n", + "ax.axhline(anomaly_fc.sel(lat=lat, lon=lon, method=\"nearest\"), color=\"k\")\n", "\n", - "ax.axhline(anomaly_fc.sel(lat=lat, lon = lon, method = 'nearest'), color = 'k')\n", + "quantiles = data.quantile([0.1, 0.9], dim=\"number\")\n", + "for q, val in zip(quantiles[\"quantile\"].values, quantiles.values):\n", + " ax.axhline(\n", + " val, color=\"k\", linestyle=\"--\", linewidth=1, label=f\"{q:.0%} quantile anom.\"\n", + " )\n", "\n", - "quantiles = data.quantile([0.1, 0.9], dim='number')\n", - "for q, val in zip(quantiles['quantile'].values, quantiles.values):\n", - " ax.axhline(val, color='k', linestyle='--', linewidth=1, label=f'{q:.0%} quantile anom.')\n", + "quantiles = years_emp.sel(lat=lat, lon=lon, method=\"nearest\").quantile(\n", + " [0.1, 0.9], dim=\"number\"\n", + ")\n", + "for q, val in zip(quantiles[\"quantile\"].values, quantiles.values):\n", + " ax.axvline(\n", + " val, color=\"red\", linestyle=\"--\", linewidth=1, label=f\"{q:.0%} quantile yr\"\n", + " )\n", "\n", - "quantiles = years_emp.sel(lat=lat, lon = lon, method = 'nearest').quantile([0.1, 0.9], dim='number')\n", - "for q, val in zip(quantiles['quantile'].values, quantiles.values):\n", - " ax.axvline(val, color='red', linestyle='--', linewidth=1, label=f'{q:.0%} quantile yr')\n", - " \n", "\n", - "ax.axvline(p_surpass.sel(lat=lat, lon = lon, method = 'nearest').where(p_surpass.sel(lat=lat, lon = lon, method = 'nearest') >= 0.5).idxmin('year'), color = 'red', label = 'crossing year')\n", + "ax.axvline(\n", + " p_surpass.sel(lat=lat, lon=lon, method=\"nearest\")\n", + " .where(p_surpass.sel(lat=lat, lon=lon, method=\"nearest\") >= 0.5)\n", + " .idxmin(\"year\"),\n", + " color=\"red\",\n", + " label=\"crossing year\",\n", + ")\n", "ax.legend()" ] }, @@ -614,7 +678,7 @@ } ], "source": [ - "ds_out['ci_width_95'].plot(cmap = 'Reds', vmin = 0)" + "ds_out[\"ci_width_95\"].plot(cmap=\"Reds\", vmin=0)" ] }, { @@ -639,7 +703,7 @@ "DATA_DIR = os.environ[\"DATA_DIR\"]\n", "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "\n", - "#Impact Regions\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(os.path.join(DATA_DIR, IMPACT_REGION_POLYGONS))\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -669,10 +733,10 @@ } ], "source": [ - "#Log forecast anomaly\n", - "fc_anomaly = isku_utils.grid_to_ir(ds2.median(dim = 'number'))\n", + "# Log forecast anomaly\n", + "fc_anomaly = isku_utils.grid_to_ir(ds2.median(dim=\"number\"))\n", "fc_anomaly_df = analysis_utils.xarray_to_gpd(fc_anomaly, _polygons)\n", - "fc_anomaly_df['value'].to_csv('2608_fc_anomaly.csv')" + "fc_anomaly_df[\"value\"].to_csv(\"2608_fc_anomaly.csv\")" ] }, { @@ -734,26 +798,26 @@ ], "source": [ "_polygons_postcard_median = analysis_utils.xarray_to_gpd(data_ir_median, _polygons)\n", - "fig, axs = plt.subplots(1,2, figsize = (12,6))\n", + "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "analysis_utils.plot_single(\n", " _polygons_postcard_median,\n", - " col='value',\n", - " cm = 'Reds',\n", - " vmin = 2020,\n", - " vmax = 2100,\n", - " cbar_label= 'Projected Decade',\n", - " ax = axs[1],\n", - " cbar_location='bottom'\n", + " col=\"value\",\n", + " cm=\"Reds\",\n", + " vmin=2020,\n", + " vmax=2100,\n", + " cbar_label=\"Projected Decade\",\n", + " ax=axs[1],\n", + " cbar_location=\"bottom\",\n", ")\n", "\n", "analysis_utils.plot_single(\n", " fc_anomaly_df,\n", - " col='value',\n", - " cm = 'bwr',\n", - " n_colors = 13,\n", - " cbar_label= 'Forecast Anomaly [C]',\n", - " ax = axs[0],\n", - " cbar_location='bottom'\n", + " col=\"value\",\n", + " cm=\"bwr\",\n", + " n_colors=13,\n", + " cbar_label=\"Forecast Anomaly [C]\",\n", + " ax=axs[0],\n", + " cbar_location=\"bottom\",\n", ")\n", "\n", "fig.savefig(\"2608_combined_anomaly_postcard\", dpi=600, bbox_inches=\"tight\")" @@ -780,27 +844,27 @@ "_polygons_postcard_mode = analysis_utils.xarray_to_gpd(data_ir_mode, _polygons)\n", "_polygons_postcard_count = analysis_utils.xarray_to_gpd(data_ir_counts, _polygons)\n", "\n", - "fig, axs = plt.subplots(1,2, figsize = (12,6))\n", + "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", "analysis_utils.plot_single(\n", " _polygons_postcard_mode,\n", - " col='value',\n", - " cm = 'Reds',\n", - " vmin = 2020,\n", - " vmax = 2100,\n", - " cbar_label= 'Projection Decade',\n", - " ax = axs[0],\n", - " cbar_location='bottom'\n", + " col=\"value\",\n", + " cm=\"Reds\",\n", + " vmin=2020,\n", + " vmax=2100,\n", + " cbar_label=\"Projection Decade\",\n", + " ax=axs[0],\n", + " cbar_location=\"bottom\",\n", ")\n", "\n", "analysis_utils.plot_single(\n", " _polygons_postcard_count,\n", - " col='value',\n", - " cm = 'Grays',\n", - " vmin = 0,\n", - " vmax = 50,\n", - " cbar_label= 'Ensemble Member Count',\n", - " ax = axs[1],\n", - " cbar_location='bottom'\n", + " col=\"value\",\n", + " cm=\"Grays\",\n", + " vmin=0,\n", + " vmax=50,\n", + " cbar_label=\"Ensemble Member Count\",\n", + " ax=axs[1],\n", + " cbar_location=\"bottom\",\n", ")\n", "\n", "fig.savefig(\"2608_combined_spread_postcard\", dpi=600, bbox_inches=\"tight\")" @@ -827,12 +891,12 @@ "_polygons_postcard_median = analysis_utils.xarray_to_gpd(data_ir_median, _polygons)\n", "ax = analysis_utils.plot_single(\n", " _polygons_postcard_median,\n", - " col='value',\n", - " cm = 'Reds',\n", - " vmin = 2020,\n", - " vmax = 2100,\n", - " cbar_label= 'Projection Year',\n", - " cbar_location='bottom'\n", + " col=\"value\",\n", + " cm=\"Reds\",\n", + " vmin=2020,\n", + " vmax=2100,\n", + " cbar_label=\"Projection Year\",\n", + " cbar_location=\"bottom\",\n", ")\n", "\n", "fig = ax.get_figure()\n", @@ -846,8 +910,10 @@ "metadata": {}, "outputs": [], "source": [ - "#df_out = _polygons_postcard[['region', 'value']].rename(columns={\"value\": \"year\"}).assign(mode_members=_polygons_postcard_count['value'])\n", - "df_out = _polygons_postcard_median[['region', 'value']].rename(columns={\"value\": \"year\"})\n", + "# df_out = _polygons_postcard[['region', 'value']].rename(columns={\"value\": \"year\"}).assign(mode_members=_polygons_postcard_count['value'])\n", + "df_out = _polygons_postcard_median[[\"region\", \"value\"]].rename(\n", + " columns={\"value\": \"year\"}\n", + ")\n", "df_out.to_csv(\"2608_postcard_median.csv\")" ] }, @@ -880,12 +946,12 @@ } ], "source": [ - "ds2 = forecast['t2m']*0+1.5\n", + "ds2 = forecast[\"t2m\"] * 0 + 1.5\n", "in_range = (ds2 >= decade_min) & (ds2 <= decade_max)\n", - "below = ds2 < decade_min.min('decade')\n", - "above = ds2 > decade_max.max('decade')\n", + "below = ds2 < decade_min.min(\"decade\")\n", + "above = ds2 > decade_max.max(\"decade\")\n", "\n", - "mode_decade = in_range.where(in_range ==True).idxmin('decade').compute()\n", + "mode_decade = in_range.where(in_range == True).idxmin(\"decade\").compute()\n", "\n", "da = xr.where(below, 2020, xr.where(above, 2110, mode_decade)).compute()\n", "data_ir = isku_utils.grid_to_ir(da)" @@ -921,12 +987,12 @@ "source": [ "ax = analysis_utils.plot_single(\n", " _polygons_postcard,\n", - " col='value',\n", - " cm = 'Reds',\n", - " vmin = 2020,\n", - " vmax = 2100,\n", - " cbar_label= 'Projection Year',\n", - " cbar_location='bottom'\n", + " col=\"value\",\n", + " cm=\"Reds\",\n", + " vmin=2020,\n", + " vmax=2100,\n", + " cbar_label=\"Projection Year\",\n", + " cbar_location=\"bottom\",\n", ")\n", "\n", "fig = ax.get_figure()\n", @@ -961,7 +1027,7 @@ } ], "source": [ - "anomaly_proj.sel(year = 2080).plot(cmap = 'RdBu_r', vmin = -5, vmax = 5)" + "anomaly_proj.sel(year=2080).plot(cmap=\"RdBu_r\", vmin=-5, vmax=5)" ] }, { @@ -982,21 +1048,21 @@ } ], "source": [ - "x = _polygons_postcard['value'][_polygons_postcard['value']>2019]\n", - "y = _polygons_postcard_count['value'][_polygons_postcard['value']>2019]\n", + "x = _polygons_postcard[\"value\"][_polygons_postcard[\"value\"] > 2019]\n", + "y = _polygons_postcard_count[\"value\"][_polygons_postcard[\"value\"] > 2019]\n", "mask = np.isfinite(x) & np.isfinite(y)\n", "\n", "fig, ax = plt.subplots()\n", - "ax.scatter(x[mask], y[mask], s=10, color='#2a78d6', alpha=0.5, edgecolors='none')\n", - "ax.set_xlabel('Mode decade')\n", - "ax.set_ylabel('Mode count')\n", - "ax.grid(color='#e1e0d9', linewidth=0.5)\n", + "ax.scatter(x[mask], y[mask], s=10, color=\"#2a78d6\", alpha=0.5, edgecolors=\"none\")\n", + "ax.set_xlabel(\"Mode decade\")\n", + "ax.set_ylabel(\"Mode count\")\n", + "ax.grid(color=\"#e1e0d9\", linewidth=0.5)\n", "ax.set_axisbelow(True)\n", - "for spine in ('top', 'right'):\n", + "for spine in (\"top\", \"right\"):\n", " ax.spines[spine].set_visible(False)\n", - "for spine in ('left', 'bottom'):\n", - " ax.spines[spine].set_color('#c3c2b7')\n", - "ax.tick_params(colors='#898781')" + "for spine in (\"left\", \"bottom\"):\n", + " ax.spines[spine].set_color(\"#c3c2b7\")\n", + "ax.tick_params(colors=\"#898781\")" ] }, { diff --git a/analysis/a08_qdm_may.ipynb b/analysis/a08_qdm_may.ipynb index e2f4fb4..dd64f89 100644 --- a/analysis/a08_qdm_may.ipynb +++ b/analysis/a08_qdm_may.ipynb @@ -31,20 +31,16 @@ ], "source": [ "import os\n", - "from pathlib import Path\n", - "from dotenv import load_dotenv, find_dotenv\n", + "from dotenv import load_dotenv\n", "import xarray as xr\n", "\n", "import geopandas as gpd\n", - "import pandas as pd\n", - "import numpy as np\n", - "import math\n", - "import matplotlib.pyplot as plt\n", "\n", "import analysis_utils\n", "import isku_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -62,7 +58,7 @@ "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", "\n", - "#Climate Data\n", + "# Climate Data\n", "TAS_FORECAST_URI = os.environ[\"POREALLAS_TAS_FORECAST_URI\"]\n", "ERA5_URI = os.environ[\"POREALLAS_ERA5_URI\"]" ] @@ -74,11 +70,14 @@ "metadata": {}, "outputs": [], "source": [ - "#Projection Effects\n", - "effect = xr.open_datatree(EFFECTS_URI, consolidated = False)\n", - "effect_qdm = xr.open_datatree(\"gs://poreallas-public-20260605/v20260731/parsed/effects.zarr\", backend_kwargs={\"storage_options\": {\"token\": \"anon\"}})\n", - "baseline_period = analysis_utils.get_baseline_period(effect, years = 30)\n", - "#Impact Regions\n", + "# Projection Effects\n", + "effect = xr.open_datatree(EFFECTS_URI, consolidated=False)\n", + "effect_qdm = xr.open_datatree(\n", + " \"gs://poreallas-public-20260605/v20260731/parsed/effects.zarr\",\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + ")\n", + "baseline_period = analysis_utils.get_baseline_period(effect, years=30)\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(IMPACT_REGION_POLYGONS)\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -86,9 +85,11 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", - "socioeconomics = socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop', 'gdppc', 'iso3']]" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", + "socioeconomics = socioeconomics.sel(year=2026)[\n", + " [\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\", \"gdppc\", \"iso3\"]\n", + "]" ] }, { @@ -98,14 +99,26 @@ "metadata": {}, "outputs": [], "source": [ - "impact_delta = analysis_utils.compute_impact(effect, socioeconomics, \n", - " chunks = {'number':-1, 'region': 'auto'},\n", - " ensemble = True, baseline_period=baseline_period, hotonly = \"coldonly\", \n", - " rate = False, age_weight = True)\n", - "impact_qdm = analysis_utils.compute_impact(effect_qdm, socioeconomics, \n", - " chunks = {'number':-1, 'region': 'auto'},\n", - " ensemble = True, baseline_period=baseline_period, hotonly = \"coldonly\", \n", - " rate = False, age_weight = True)" + "impact_delta = analysis_utils.compute_impact(\n", + " effect,\n", + " socioeconomics,\n", + " chunks={\"number\": -1, \"region\": \"auto\"},\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=\"coldonly\",\n", + " rate=False,\n", + " age_weight=True,\n", + ")\n", + "impact_qdm = analysis_utils.compute_impact(\n", + " effect_qdm,\n", + " socioeconomics,\n", + " chunks={\"number\": -1, \"region\": \"auto\"},\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=\"coldonly\",\n", + " rate=False,\n", + " age_weight=True,\n", + ")" ] }, { @@ -144,12 +157,31 @@ } ], "source": [ - "#QDM\n", - "era5_qdm = xr.open_dataset(\"gs://poreallas-public-20260605/v20260731/parsed/era5_adj.zarr\", chunks = {'time':-1}, backend_kwargs={\"storage_options\": {\"token\": \"anon\"}}).groupby(\"time.month\").mean().sel(month = [5, 6, 7, 8, 9, 10])\n", - "forecast_qdm = xr.open_dataset(\"gs://poreallas-public-20260605/v20260731/parsed/forecast_adj.zarr\", chunks = {'time':-1}, backend_kwargs={\"storage_options\": {\"token\": \"anon\"}}).mean(dim = 'number').groupby(\"time.month\").mean().sel(month = [5, 6, 7, 8, 9, 10])\n", + "# QDM\n", + "era5_qdm = (\n", + " xr.open_dataset(\n", + " \"gs://poreallas-public-20260605/v20260731/parsed/era5_adj.zarr\",\n", + " chunks={\"time\": -1},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[5, 6, 7, 8, 9, 10])\n", + ")\n", + "forecast_qdm = (\n", + " xr.open_dataset(\n", + " \"gs://poreallas-public-20260605/v20260731/parsed/forecast_adj.zarr\",\n", + " chunks={\"time\": -1},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )\n", + " .mean(dim=\"number\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[5, 6, 7, 8, 9, 10])\n", + ")\n", "\n", - "forecast_qdm = isku_utils.grid_to_ir(forecast_qdm, savefile = None)\n", - "era5_qdm = isku_utils.grid_to_ir(era5_qdm, savefile = None)" + "forecast_qdm = isku_utils.grid_to_ir(forecast_qdm, savefile=None)\n", + "era5_qdm = isku_utils.grid_to_ir(era5_qdm, savefile=None)" ] }, { @@ -160,8 +192,18 @@ "outputs": [], "source": [ "## Delta Shift\n", - "forecast_gmfd = xr.open_zarr(os.path.join(DATA_DIR, \"daily_forecast_GMFD_tas.zarr\")).mean(dim = 'number').groupby(\"time.month\").mean()\n", - "era5_gmfd = xr.open_zarr(os.path.join(DATA_DIR, \"daily_reanalysis_GMFD_tas.zarr\")).groupby('time.month').mean().sel(month = [5, 6, 7, 8, 9, 10])" + "forecast_gmfd = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"daily_forecast_GMFD_tas.zarr\"))\n", + " .mean(dim=\"number\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + ")\n", + "era5_gmfd = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"daily_reanalysis_GMFD_tas.zarr\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[5, 6, 7, 8, 9, 10])\n", + ")" ] }, { @@ -192,12 +234,23 @@ } ], "source": [ - "#Raw\n", - "forecast_raw = xr.open_zarr(os.path.join(DATA_DIR, \"tas_forecast.zarr\")).mean(dim = 'number').swap_dims({'forecast_period': 'valid_time'}).groupby(\"valid_time.month\").mean()\n", - "era5_raw = xr.open_zarr(os.path.join(DATA_DIR, \"tas_ERA5.zarr\")).groupby('time.month').mean().sel(month = [5, 6, 7, 8, 9, 10])\n", + "# Raw\n", + "forecast_raw = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"tas_forecast.zarr\"))\n", + " .mean(dim=\"number\")\n", + " .swap_dims({\"forecast_period\": \"valid_time\"})\n", + " .groupby(\"valid_time.month\")\n", + " .mean()\n", + ")\n", + "era5_raw = (\n", + " xr.open_zarr(os.path.join(DATA_DIR, \"tas_ERA5.zarr\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[5, 6, 7, 8, 9, 10])\n", + ")\n", "\n", - "forecast_raw = isku_utils.grid_to_ir(forecast_raw, savefile = None)\n", - "era5_raw = isku_utils.grid_to_ir(era5_raw, savefile = None)" + "forecast_raw = isku_utils.grid_to_ir(forecast_raw, savefile=None)\n", + "era5_raw = isku_utils.grid_to_ir(era5_raw, savefile=None)" ] }, { @@ -218,19 +271,20 @@ } ], "source": [ - "diff = forecast_qdm-forecast_gmfd\n", + "diff = forecast_qdm - forecast_gmfd\n", "\n", - "_polygons_delta = analysis_utils.xarray_to_gpd(diff['tas'], _polygons)\n", + "_polygons_delta = analysis_utils.xarray_to_gpd(diff[\"tas\"], _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='tas',\n", - " cm = 'bwr',\n", - " vmin = -5, vmax = 5,\n", - " sup_title=f\"\",\n", - " cbar_label= \"anomaly_forecast_qdm_delta\",\n", + " col=\"tas\",\n", + " cm=\"bwr\",\n", + " vmin=-5,\n", + " vmax=5,\n", + " sup_title=\"\",\n", + " cbar_label=\"anomaly_forecast_qdm_delta\",\n", ")\n", "\n", - "fig.savefig(f\"anomaly_forecast_qdm_delta.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\"anomaly_forecast_qdm_delta.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -251,19 +305,20 @@ } ], "source": [ - "diff = forecast_qdm-forecast_raw\n", + "diff = forecast_qdm - forecast_raw\n", "\n", - "_polygons_delta = analysis_utils.xarray_to_gpd(diff['tas'], _polygons)\n", + "_polygons_delta = analysis_utils.xarray_to_gpd(diff[\"tas\"], _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='tas',\n", - " cm = 'bwr',\n", - " vmin = -5, vmax = 5,\n", - " sup_title=f\"\",\n", - " cbar_label= \"anomaly_forecast_qdm_raw\",\n", + " col=\"tas\",\n", + " cm=\"bwr\",\n", + " vmin=-5,\n", + " vmax=5,\n", + " sup_title=\"\",\n", + " cbar_label=\"anomaly_forecast_qdm_raw\",\n", ")\n", "\n", - "fig.savefig(f\"anomaly_forecast_qdm_raw.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\"anomaly_forecast_qdm_raw.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -292,19 +347,20 @@ } ], "source": [ - "diff = era5_qdm-era5_gmfd\n", + "diff = era5_qdm - era5_gmfd\n", "\n", - "_polygons_delta = analysis_utils.xarray_to_gpd(diff['tas'], _polygons)\n", + "_polygons_delta = analysis_utils.xarray_to_gpd(diff[\"tas\"], _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='tas',\n", - " cm = 'bwr',\n", - " vmin = -5, vmax = 5,\n", - " sup_title=f\"\",\n", - " cbar_label= \"anomaly_era5_qdm_delta\",\n", + " col=\"tas\",\n", + " cm=\"bwr\",\n", + " vmin=-5,\n", + " vmax=5,\n", + " sup_title=\"\",\n", + " cbar_label=\"anomaly_era5_qdm_delta\",\n", ")\n", "\n", - "fig.savefig(f\"anomaly_era5_qdm_delta.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\"anomaly_era5_qdm_delta.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -325,19 +381,20 @@ } ], "source": [ - "diff = era5_qdm-era5_raw\n", + "diff = era5_qdm - era5_raw\n", "\n", - "_polygons_delta = analysis_utils.xarray_to_gpd(diff['tas'], _polygons)\n", + "_polygons_delta = analysis_utils.xarray_to_gpd(diff[\"tas\"], _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='tas',\n", - " cm = 'bwr',\n", - " vmin = -5, vmax = 5,\n", - " sup_title=f\"\",\n", - " cbar_label= \"anomaly_era5_qdm_raw\",\n", + " col=\"tas\",\n", + " cm=\"bwr\",\n", + " vmin=-5,\n", + " vmax=5,\n", + " sup_title=\"\",\n", + " cbar_label=\"anomaly_era5_qdm_raw\",\n", ")\n", "\n", - "fig.savefig(f\"anomaly_era5_qdm_raw.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\"anomaly_era5_qdm_raw.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -355,19 +412,20 @@ "metadata": {}, "outputs": [], "source": [ - "diff = forecast_qdm-era5_qdm\n", + "diff = forecast_qdm - era5_qdm\n", "\n", - "_polygons_delta = analysis_utils.xarray_to_gpd(diff['tas'], _polygons)\n", + "_polygons_delta = analysis_utils.xarray_to_gpd(diff[\"tas\"], _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='tas',\n", - " cm = 'bwr',\n", - " vmin = -5, vmax = 5,\n", - " sup_title=f\"\",\n", - " cbar_label= \"anomaly_forecast_era5_qdm\",\n", + " col=\"tas\",\n", + " cm=\"bwr\",\n", + " vmin=-5,\n", + " vmax=5,\n", + " sup_title=\"\",\n", + " cbar_label=\"anomaly_forecast_era5_qdm\",\n", ")\n", "\n", - "fig.savefig(f\"anomaly_forecast_era5_qdm.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\"anomaly_forecast_era5_qdm.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -396,22 +454,35 @@ } ], "source": [ - "effect_delta_temp = effect['/baseline']['effect'].sel(age_cohort = 'age65plus').groupby('time.month').mean()\n", - "effect_qdm_temp = effect_qdm['/baseline']['effect'].sel(age_cohort = 'age65plus').groupby('time.month').mean()\n", + "effect_delta_temp = (\n", + " effect[\"/baseline\"][\"effect\"]\n", + " .sel(age_cohort=\"age65plus\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + ")\n", + "effect_qdm_temp = (\n", + " effect_qdm[\"/baseline\"][\"effect\"]\n", + " .sel(age_cohort=\"age65plus\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + ")\n", "\n", - "diff = effect_qdm_temp-effect_delta_temp\n", + "diff = effect_qdm_temp - effect_delta_temp\n", "\n", "_polygons_delta = analysis_utils.xarray_to_gpd(diff, _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='effect',\n", - " cm = 'bwr',\n", - " vmax = 20, vmin = -20,\n", - " sup_title=f\"\",\n", - " cbar_label= \"diff_mortality_rate_effect\",\n", + " col=\"effect\",\n", + " cm=\"bwr\",\n", + " vmax=20,\n", + " vmin=-20,\n", + " sup_title=\"\",\n", + " cbar_label=\"diff_mortality_rate_effect\",\n", ")\n", "\n", - "fig.savefig(f\"diff_mortality_rate_65plus_effect_baseline.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\n", + " \"diff_mortality_rate_65plus_effect_baseline.png\", dpi=600, bbox_inches=\"tight\"\n", + ")" ] }, { @@ -432,22 +503,37 @@ } ], "source": [ - "effect_delta_temp = effect['/forecast']['effect'].sel(age_cohort = 'age65plus').groupby('time.month').mean().mean(dim = 'number')\n", - "effect_qdm_temp = effect_qdm['/forecast']['effect'].sel(age_cohort = 'age65plus').groupby('time.month').mean().mean(dim = 'number')\n", + "effect_delta_temp = (\n", + " effect[\"/forecast\"][\"effect\"]\n", + " .sel(age_cohort=\"age65plus\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .mean(dim=\"number\")\n", + ")\n", + "effect_qdm_temp = (\n", + " effect_qdm[\"/forecast\"][\"effect\"]\n", + " .sel(age_cohort=\"age65plus\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .mean(dim=\"number\")\n", + ")\n", "\n", - "diff = effect_qdm_temp-effect_delta_temp\n", + "diff = effect_qdm_temp - effect_delta_temp\n", "\n", "_polygons_delta = analysis_utils.xarray_to_gpd(diff, _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='effect',\n", - " cm = 'bwr',\n", - " vmax = 20, vmin = -20,\n", - " sup_title=f\"\",\n", - " cbar_label= \"diff_mortality_rate_effect\",\n", + " col=\"effect\",\n", + " cm=\"bwr\",\n", + " vmax=20,\n", + " vmin=-20,\n", + " sup_title=\"\",\n", + " cbar_label=\"diff_mortality_rate_effect\",\n", ")\n", "\n", - "fig.savefig(f\"diff_mortality_rate_65plus_effect_forecast.png\", dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\n", + " \"diff_mortality_rate_65plus_effect_forecast.png\", dpi=600, bbox_inches=\"tight\"\n", + ")" ] }, { @@ -484,14 +570,17 @@ } ], "source": [ - "_polygons_delta = analysis_utils.xarray_to_gpd(impact_delta.mean(dim = 'number'), _polygons)\n", + "_polygons_delta = analysis_utils.xarray_to_gpd(\n", + " impact_delta.mean(dim=\"number\"), _polygons\n", + ")\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_delta,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " vmax = 7, vmin = -7,\n", - " sup_title=f\"\",\n", - " cbar_label= \"delta_mortality_rate\",\n", + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " vmax=7,\n", + " vmin=-7,\n", + " sup_title=\"\",\n", + " cbar_label=\"delta_mortality_rate\",\n", ")\n", "\n", "fig.savefig(\"delta_mortality_rate.png\", dpi=600, bbox_inches=\"tight\")" @@ -523,14 +612,15 @@ } ], "source": [ - "_polygons_qdm = analysis_utils.xarray_to_gpd(impact_qdm.mean(dim = 'number'), _polygons)\n", + "_polygons_qdm = analysis_utils.xarray_to_gpd(impact_qdm.mean(dim=\"number\"), _polygons)\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_qdm,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " vmax = 7, vmin = -7,\n", - " sup_title=f\"\",\n", - " cbar_label= \"qdm_mortality_rate\",\n", + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " vmax=7,\n", + " vmin=-7,\n", + " sup_title=\"\",\n", + " cbar_label=\"qdm_mortality_rate\",\n", ")\n", "\n", "fig.savefig(\"qdm_mortality_rate.png\", dpi=600, bbox_inches=\"tight\")" @@ -554,16 +644,17 @@ } ], "source": [ - "diff = impact_qdm-impact_delta\n", - "_polygons_impact = analysis_utils.xarray_to_gpd(diff.mean(dim = 'number'), _polygons)\n", + "diff = impact_qdm - impact_delta\n", + "_polygons_impact = analysis_utils.xarray_to_gpd(diff.mean(dim=\"number\"), _polygons)\n", "\n", "fig = analysis_utils.plot_monthly(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " vmax = 7, vmin = -7,\n", - " sup_title=f\"\",\n", - " cbar_label= \"diff_mortality_rate (qdm-delta)\",\n", + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " vmax=7,\n", + " vmin=-7,\n", + " sup_title=\"\",\n", + " cbar_label=\"diff_mortality_rate (qdm-delta)\",\n", ")\n", "\n", "fig.savefig(\"diff_mortality_rate_qdmSUBdelta.png\", dpi=600, bbox_inches=\"tight\")" diff --git a/analysis/a10_Effect_Diagnostic_Panels.ipynb b/analysis/a10_Effect_Diagnostic_Panels.ipynb index a082401..01b3d29 100644 --- a/analysis/a10_Effect_Diagnostic_Panels.ipynb +++ b/analysis/a10_Effect_Diagnostic_Panels.ipynb @@ -42,21 +42,19 @@ "\n", "import geopandas as gpd\n", "import pandas as pd\n", - "import marimo as mo\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import seaborn as sns\n", "import xarray as xr\n", "import math\n", "import matplotlib.colors as mcolors\n", "\n", - "import cartopy.crs as ccrs\n", "\n", "import analysis_utils\n", "import isku_utils\n", "import diagnostic_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -82,7 +80,7 @@ "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", "\n", - "#Climate Data\n", + "# Climate Data\n", "TAS_FORECAST_URI = os.environ[\"POREALLAS_TAS_FORECAST_URI\"]\n", "ERA5_URI = os.environ[\"POREALLAS_ERA5_URI\"]" ] @@ -94,9 +92,9 @@ "metadata": {}, "outputs": [], "source": [ - "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated = False)\n", - "baseline_period = analysis_utils.get_baseline_period(effect, years = 30)\n", - "#Impact Regions\n", + "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated=False)\n", + "baseline_period = analysis_utils.get_baseline_period(effect, years=30)\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(IMPACT_REGION_POLYGONS)\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -104,9 +102,11 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", - "socioeconomics = socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop', 'gdppc', 'iso3']]" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", + "socioeconomics = socioeconomics.sel(year=2026)[\n", + " [\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\", \"gdppc\", \"iso3\"]\n", + "]" ] }, { @@ -155,12 +155,18 @@ "source": [ "# Climate Data by Impact Regoion\n", "\n", - "forecast = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/2608_forecast_adj_corrected.zarr\")\n", - "forecast_ir = isku_utils.grid_to_ir(forecast, savefile = None)\n", - "reanalysis = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/era5_adj_corrected.zarr\")\n", - "reanalysis_ir = isku_utils.grid_to_ir(reanalysis, savefile = None)\n", - "#Intermediate Betas/MMT from Project.Py\n", - "betas_mmt = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/2608_beta_qdm.zarr\")\n", + "forecast = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/2608_forecast_adj_corrected.zarr\"\n", + ")\n", + "forecast_ir = isku_utils.grid_to_ir(forecast, savefile=None)\n", + "reanalysis = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/era5_adj_corrected.zarr\"\n", + ")\n", + "reanalysis_ir = isku_utils.grid_to_ir(reanalysis, savefile=None)\n", + "# Intermediate Betas/MMT from Project.Py\n", + "betas_mmt = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/2608_beta_qdm.zarr\"\n", + ")\n", "\n", "projected = effect" ] @@ -183,18 +189,15 @@ } ], "source": [ - "import matplotlib.patheffects as pe\n", "from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec\n", "\n", - "import cartopy.crs as ccrs\n", - "from scipy.interpolate import interp1d\n", "from shapely import box\n", "\n", "### Analysis Options ###\n", - "rate = True # Mortality Rate vs Total Deaths\n", - "hotonly = False # Hot-Only Deaths\n", - "age_weight = False #Proportional Age-weighting vs single age group\n", - "density_hist = True #Controls daily data histograms (not computation normalization)\n", + "rate = True # Mortality Rate vs Total Deaths\n", + "hotonly = False # Hot-Only Deaths\n", + "age_weight = False # Proportional Age-weighting vs single age group\n", + "density_hist = True # Controls daily data histograms (not computation normalization)\n", "\n", "### Region Selection ###\n", "region_filters = [\n", @@ -213,13 +216,18 @@ "]\n", "\n", "for region_filter in [\"ITA.8.39\"]:\n", - "\n", " # Compute Impact\n", - " regional_sum = analysis_utils.compute_impact(projected, socioeconomics, baseline_period= baseline_period, \n", - " hotonly = hotonly, rate = rate, age_weight = age_weight,\n", - " chunks = {'number':-1, 'region': 'auto'})\n", - " \n", - " #Get Months\n", + " regional_sum = analysis_utils.compute_impact(\n", + " projected,\n", + " socioeconomics,\n", + " baseline_period=baseline_period,\n", + " hotonly=hotonly,\n", + " rate=rate,\n", + " age_weight=age_weight,\n", + " chunks={\"number\": -1, \"region\": \"auto\"},\n", + " )\n", + "\n", + " # Get Months\n", " months = regional_sum.month.values[1:]\n", " base_month = months[0]\n", " lead_by_month = {month: (month - base_month + 1) for month in months}\n", @@ -232,25 +240,45 @@ " _country_proj = _country.to_crs(\"ESRI:54030\")\n", " bminx, bminy, bmaxx, bmaxy = _country_proj.total_bounds\n", " pad_proj = 50000\n", - " \n", + "\n", " # BBOX Impact Regions\n", - " _polygons_local = _polygons.to_crs(\"ESRI:54030\").clip(box(bminx-pad_proj, bminy-pad_proj, bmaxx+pad_proj, bmaxy+pad_proj)).reset_index()\n", + " _polygons_local = (\n", + " _polygons.to_crs(\"ESRI:54030\")\n", + " .clip(\n", + " box(bminx - pad_proj, bminy - pad_proj, bmaxx + pad_proj, bmaxy + pad_proj)\n", + " )\n", + " .reset_index()\n", + " )\n", "\n", " ### Load Climate Data ###\n", - " forecast_local = forecast_ir['tas'].sel(region = [region_filter])-273.15\n", - " reanalysis_local = reanalysis_ir['tas'].sel(region = [region_filter])-273.15\n", - " tas_anomaly = forecast_local.groupby(\"time.month\").mean().mean(dim = 'number')-reanalysis_local.groupby(\"time.month\").mean()\n", - "\n", - " combined = np.concat([forecast_local.values.flatten(), reanalysis_local.values.flatten()])\n", + " forecast_local = forecast_ir[\"tas\"].sel(region=[region_filter]) - 273.15\n", + " reanalysis_local = reanalysis_ir[\"tas\"].sel(region=[region_filter]) - 273.15\n", + " tas_anomaly = (\n", + " forecast_local.groupby(\"time.month\").mean().mean(dim=\"number\")\n", + " - reanalysis_local.groupby(\"time.month\").mean()\n", + " )\n", + "\n", + " combined = np.concat(\n", + " [forecast_local.values.flatten(), reanalysis_local.values.flatten()]\n", + " )\n", " xmin, xmax = np.percentile(combined, [1, 100])\n", "\n", " # Precompute Impact\n", - " regional_sum_roi = regional_sum.sel(region = region_filter)\n", + " regional_sum_roi = regional_sum.sel(region=region_filter)\n", " absmax_impact = np.abs(regional_sum_roi).quantile(0.95).values\n", - " norm_impact = mcolors.TwoSlopeNorm(vmin=-absmax_impact, vcenter=0, vmax=absmax_impact)\n", + " norm_impact = mcolors.TwoSlopeNorm(\n", + " vmin=-absmax_impact, vcenter=0, vmax=absmax_impact\n", + " )\n", "\n", " # Compute Cumulative Effect\n", - " cdf_data, max_cdf = diagnostic_utils.compute_cumulative_effect(forecast_local, reanalysis_local, region_filter, months, monthly = True, hotonly = hotonly)\n", + " cdf_data, max_cdf = diagnostic_utils.compute_cumulative_effect(\n", + " forecast_local,\n", + " reanalysis_local,\n", + " region_filter,\n", + " months,\n", + " monthly=True,\n", + " hotonly=hotonly,\n", + " )\n", "\n", " # Initialize Figure\n", " ncols = 3\n", @@ -266,33 +294,38 @@ " hist_axes = []\n", "\n", " for pair_row in range(n_pairs):\n", - " inner_gs = GridSpecFromSubplotSpec(3, ncols, subplot_spec=gs[pair_row],\n", - " hspace=0.15, wspace=0.25,\n", - " height_ratios=[2, 2, 1.2])\n", + " inner_gs = GridSpecFromSubplotSpec(\n", + " 3,\n", + " ncols,\n", + " subplot_spec=gs[pair_row],\n", + " hspace=0.15,\n", + " wspace=0.25,\n", + " height_ratios=[2, 2, 1.2],\n", + " )\n", " for col in range(ncols):\n", " i = pair_row * ncols + col\n", " if i >= n_months:\n", " break\n", " ax_top = fig.add_subplot(inner_gs[0, col])\n", " ax_bot = fig.add_subplot(inner_gs[1, col])\n", - " hist_gs = GridSpecFromSubplotSpec(1, 3, subplot_spec=inner_gs[2, col], width_ratios=[0.2, 0.6, 0.2])\n", + " hist_gs = GridSpecFromSubplotSpec(\n", + " 1, 3, subplot_spec=inner_gs[2, col], width_ratios=[0.2, 0.6, 0.2]\n", + " )\n", " ax_hist = fig.add_subplot(hist_gs[0, 1])\n", " top_axes.append(ax_top)\n", " bot_axes.append(ax_bot)\n", " hist_axes.append(ax_hist)\n", "\n", - " \n", - "\n", " ### Plot ###\n", " for i, month in enumerate(months):\n", " ax_top = top_axes[i]\n", " ax_bot = bot_axes[i]\n", " ax_hist = hist_axes[i]\n", - " month_label = pd.Timestamp(2026, month, 1).strftime('%B')\n", + " month_label = pd.Timestamp(2026, month, 1).strftime(\"%B\")\n", "\n", " # Top: temperature\n", " _polygons_data_tas = _polygons_local.merge(\n", - " tas_anomaly.sel(month = month).to_dataframe().reset_index(),\n", + " tas_anomaly.sel(month=month).to_dataframe().reset_index(),\n", " on=\"region\",\n", " )\n", "\n", @@ -302,16 +335,19 @@ " legend=False,\n", " ax=ax_top,\n", " cmap=\"RdBu_r\",\n", - " vmin=-5, vmax=5,\n", - " edgecolor='black',\n", + " vmin=-5,\n", + " vmax=5,\n", + " edgecolor=\"black\",\n", " linewidth=0.03,\n", " )\n", - " _plot_data_forecast[_plot_data_forecast[\"region\"] == region_filter].boundary.plot(ax=ax_top, color=\"black\", linewidth=1.5)\n", + " _plot_data_forecast[\n", + " _plot_data_forecast[\"region\"] == region_filter\n", + " ].boundary.plot(ax=ax_top, color=\"black\", linewidth=1.5)\n", " ax_top.set_title(month_label)\n", " ax_top.set_axis_off()\n", "\n", " # Middle: impact\n", - " regional_sum_month = regional_sum.sel(month = month)\n", + " regional_sum_month = regional_sum.sel(month=month)\n", " _polygons_data = _polygons_local.merge(\n", " regional_sum_month.to_dataframe().reset_index(),\n", " on=\"region\",\n", @@ -323,25 +359,39 @@ " ax=ax_bot,\n", " cmap=\"seismic\",\n", " norm=norm_impact,\n", - " edgecolor='black',\n", + " edgecolor=\"black\",\n", " linewidth=0.03,\n", " )\n", - " _plot_data[_plot_data[\"region\"] == region_filter].boundary.plot(ax=ax_bot, color=\"black\", linewidth=1.5)\n", + " _plot_data[_plot_data[\"region\"] == region_filter].boundary.plot(\n", + " ax=ax_bot, color=\"black\", linewidth=1.5\n", + " )\n", " ax_bot.set_axis_off()\n", "\n", " # Bottom: histogram\n", - " forecast_local_month = forecast_local.sel(region = region_filter, time=forecast_local.time.dt.month == month)\n", - " reanalysis_local_month = reanalysis_local.sel(region = region_filter, time=reanalysis_local.time.dt.month == month)\n", + " forecast_local_month = forecast_local.sel(\n", + " region=region_filter, time=forecast_local.time.dt.month == month\n", + " )\n", + " reanalysis_local_month = reanalysis_local.sel(\n", + " region=region_filter, time=reanalysis_local.time.dt.month == month\n", + " )\n", "\n", - " ref_vals = betas_mmt['mmt'].sel(region=region_filter).sel(age_cohort='age65plus')\n", - " betas = betas_mmt['beta_hotonly'] if hotonly else betas_mmt['beta']\n", - " da_temp_bins = betas.sel(region=region_filter).sel(age_cohort='age65plus')\n", - " da_temp_bins['tas_bin'] = da_temp_bins['tas_bin'] \n", + " ref_vals = (\n", + " betas_mmt[\"mmt\"].sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + " )\n", + " betas = betas_mmt[\"beta_hotonly\"] if hotonly else betas_mmt[\"beta\"]\n", + " da_temp_bins = betas.sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + " da_temp_bins[\"tas_bin\"] = da_temp_bins[\"tas_bin\"]\n", "\n", " # Density histograms on hidden twin axis\n", - " bins = da_temp_bins.coords['tas_bin']\n", - " ax_hist.plot(da_temp_bins.coords['tas_bin'], da_temp_bins.values, color='red', linewidth=1.5, label=\"Response curve\")\n", - " ax_hist.axvline(ref_vals.values, color='red', linewidth=0.75, label=\"MMT\")\n", + " bins = da_temp_bins.coords[\"tas_bin\"]\n", + " ax_hist.plot(\n", + " da_temp_bins.coords[\"tas_bin\"],\n", + " da_temp_bins.values,\n", + " color=\"red\",\n", + " linewidth=1.5,\n", + " label=\"Response curve\",\n", + " )\n", + " ax_hist.axvline(ref_vals.values, color=\"red\", linewidth=0.75, label=\"MMT\")\n", " ax_hist.set_ylabel(\"Response 65+\")\n", " ax_hist.set_ylim(-10, 30)\n", " ax_hist.set_xlabel(\"Tas (C)\")\n", @@ -349,10 +399,34 @@ "\n", " # Density histograms on hidden twin axis\n", " ax_hist2 = ax_hist.twinx()\n", - " ax_hist2.hist(reanalysis_local_month.values.flatten(), bins=bins, alpha=0.5, label=\"ERA5 (1995–2025)\", density=density_hist, color='black')\n", - " ax_hist2.hist(forecast_local_month.values.flatten(), bins=bins, alpha=0.5, label=\"Forecast (51-mem.)\", density=density_hist, color='tab:orange')\n", - " ax_hist2.axvline(forecast_local_month.values.flatten().mean(), color='darkorange', linewidth=1, linestyle='--')\n", - " ax_hist2.axvline(reanalysis_local_month.values.flatten().mean(), color='black', linewidth=1, linestyle='--')\n", + " ax_hist2.hist(\n", + " reanalysis_local_month.values.flatten(),\n", + " bins=bins,\n", + " alpha=0.5,\n", + " label=\"ERA5 (1995–2025)\",\n", + " density=density_hist,\n", + " color=\"black\",\n", + " )\n", + " ax_hist2.hist(\n", + " forecast_local_month.values.flatten(),\n", + " bins=bins,\n", + " alpha=0.5,\n", + " label=\"Forecast (51-mem.)\",\n", + " density=density_hist,\n", + " color=\"tab:orange\",\n", + " )\n", + " ax_hist2.axvline(\n", + " forecast_local_month.values.flatten().mean(),\n", + " color=\"darkorange\",\n", + " linewidth=1,\n", + " linestyle=\"--\",\n", + " )\n", + " ax_hist2.axvline(\n", + " reanalysis_local_month.values.flatten().mean(),\n", + " color=\"black\",\n", + " linewidth=1,\n", + " linestyle=\"--\",\n", + " )\n", " ax_hist2.set_yticks([])\n", " ax_hist2.spines[\"right\"].set_visible(False)\n", " ax_hist2.set_xlim(xmin, xmax)\n", @@ -364,24 +438,37 @@ " ax3.spines[\"right\"].set_position((\"axes\", 1.0))\n", " ax3.set_ylabel(\"Cum. Deaths\\n/100,000\")\n", " ax3.set_ylim(0, max_cdf)\n", - " ax3.plot(fc_centers, fc_cdf, color='tab:orange', label=\"Forecast weighted CDF\")\n", - " ax3.plot(era5_centers, era5_cdf, color='black', label=\"ERA5 weighted CDF\")\n", + " ax3.plot(fc_centers, fc_cdf, color=\"tab:orange\", label=\"Forecast weighted CDF\")\n", + " ax3.plot(era5_centers, era5_cdf, color=\"black\", label=\"ERA5 weighted CDF\")\n", "\n", " if i == 0:\n", - " handles = (ax_hist.get_legend_handles_labels()[0]\n", - " + ax_hist2.get_legend_handles_labels()[0]\n", - " + ax3.get_legend_handles_labels()[0])\n", - " labels = (ax_hist.get_legend_handles_labels()[1]\n", - " + ax_hist2.get_legend_handles_labels()[1]\n", - " + ax3.get_legend_handles_labels()[1])\n", - " ax_hist.legend(handles, labels, fontsize=7, loc='upper center',\n", - " bbox_to_anchor=(0.5, -0.75), ncol=3, frameon=False)\n", + " handles = (\n", + " ax_hist.get_legend_handles_labels()[0]\n", + " + ax_hist2.get_legend_handles_labels()[0]\n", + " + ax3.get_legend_handles_labels()[0]\n", + " )\n", + " labels = (\n", + " ax_hist.get_legend_handles_labels()[1]\n", + " + ax_hist2.get_legend_handles_labels()[1]\n", + " + ax3.get_legend_handles_labels()[1]\n", + " )\n", + " ax_hist.legend(\n", + " handles,\n", + " labels,\n", + " fontsize=7,\n", + " loc=\"upper center\",\n", + " bbox_to_anchor=(0.5, -0.75),\n", + " ncol=3,\n", + " frameon=False,\n", + " )\n", "\n", " if col == ncols - 1 and pair_row == n_pairs - 1:\n", " ax_map = fig.add_subplot(inner_gs[:, col])\n", " pad_proj = 5000000\n", - " _polygons.to_crs(\"ESRI:54030\").plot(ax=ax_map, facecolor='lightgray', edgecolor='gray', linewidth=0.5)\n", - " _country_proj.plot(ax=ax_map, facecolor='gray', edgecolor='black', linewidth=1)\n", + " _polygons.to_crs(\"ESRI:54030\").plot(\n", + " ax=ax_map, facecolor=\"lightgray\", edgecolor=\"gray\", linewidth=0.5\n", + " )\n", + " _country_proj.plot(ax=ax_map, facecolor=\"gray\", edgecolor=\"black\", linewidth=1)\n", " ax_map.set_xlim(bminx - pad_proj, bmaxx + pad_proj)\n", " ax_map.set_ylim(bminy - pad_proj, bmaxy + pad_proj)\n", " ax_map.set_xticks([])\n", @@ -395,33 +482,41 @@ " fig.canvas.draw()\n", "\n", " pos_top = ax_oct_top.get_position()\n", - " cax_temp = fig.add_axes([pos_top.x1 + 0.01, pos_top.y0+0.05, 0.01, pos_top.height])\n", - " sm_temp = plt.cm.ScalarMappable(cmap='RdBu_r', norm=plt.Normalize(vmin=-5, vmax=5))\n", + " cax_temp = fig.add_axes(\n", + " [pos_top.x1 + 0.01, pos_top.y0 + 0.05, 0.01, pos_top.height]\n", + " )\n", + " sm_temp = plt.cm.ScalarMappable(cmap=\"RdBu_r\", norm=plt.Normalize(vmin=-5, vmax=5))\n", " fig.colorbar(sm_temp, cax=cax_temp, label=\"ΔTas (K)\\nBaseline: 1995-2025\")\n", "\n", " pos_bot = ax_oct_bot.get_position()\n", - " cax_impact = fig.add_axes([pos_bot.x1 + 0.01, pos_bot.y0+.05, 0.01, pos_bot.height])\n", + " cax_impact = fig.add_axes(\n", + " [pos_bot.x1 + 0.01, pos_bot.y0 + 0.05, 0.01, pos_bot.height]\n", + " )\n", " sm_impact = plt.cm.ScalarMappable(cmap=\"seismic\", norm=norm_impact)\n", " label01 = \"Monthly Impact\"\n", " units = \"[deaths/100k]\" if rate else \"[num. deaths]\"\n", - " label = f\"{label01}:\\n{projected['/forecast']['effect'].attrs['long_name']}\\n{units}\"\n", + " label = (\n", + " f\"{label01}:\\n{projected['/forecast']['effect'].attrs['long_name']}\\n{units}\"\n", + " )\n", " fig.colorbar(sm_impact, cax=cax_impact, label=label)\n", "\n", - "\n", " fig.subplots_adjust(top=0.93)\n", "\n", " title_region = f\" {region_filter}\" if region_filter else \"\"\n", " ag_label = \"Age-Weighted\" if age_weight else \"Age Cohort: 65+\"\n", - " label = 'Death Rate' if rate else \"Deaths\"\n", - " label_h = 'Hotonly' if hotonly else 'Net'\n", + " label = \"Death Rate\" if rate else \"Deaths\"\n", + " label_h = \"Hotonly\" if hotonly else \"Net\"\n", "\n", " fig.suptitle(f\"{ag_label} {label_h} {label}:{title_region}\", fontsize=14)\n", "\n", - "\n", - " label_f = 'death_rate' if rate else \"deaths\"\n", + " label_f = \"death_rate\" if rate else \"deaths\"\n", " ag_label = \"agwtdd\" if age_weight else \"65plus\"\n", "\n", - " fig.savefig(f'{label_h}_{label_f}_{ag_label}{f'_{region_filter}' if region_filter else ''}.png', dpi=600, bbox_inches=\"tight\")" + " fig.savefig(\n", + " f\"{label_h}_{label_f}_{ag_label}{f'_{region_filter}' if region_filter else ''}.png\",\n", + " dpi=600,\n", + " bbox_inches=\"tight\",\n", + " )" ] }, { diff --git a/analysis/a11_mmt_covariates.ipynb b/analysis/a11_mmt_covariates.ipynb index 5cac360..7339991 100644 --- a/analysis/a11_mmt_covariates.ipynb +++ b/analysis/a11_mmt_covariates.ipynb @@ -31,20 +31,16 @@ ], "source": [ "import os\n", - "from pathlib import Path\n", - "from dotenv import load_dotenv, find_dotenv\n", + "from dotenv import load_dotenv\n", "import xarray as xr\n", "\n", "import geopandas as gpd\n", - "import pandas as pd\n", - "import numpy as np\n", - "import math\n", - "import matplotlib.pyplot as plt\n", "\n", "import analysis_utils\n", "import isku_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -59,8 +55,7 @@ "load_dotenv()\n", "DATA_DIR = os.environ[\"DATA_DIR\"]\n", "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", - "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", - "\n" + "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]" ] }, { @@ -70,7 +65,7 @@ "metadata": {}, "outputs": [], "source": [ - "#Impact Regions\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(IMPACT_REGION_POLYGONS)\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -78,8 +73,8 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)" ] }, { @@ -1292,18 +1287,22 @@ } ], "source": [ - "_polygons_mmt = analysis_utils.xarray_to_gpd(betas['mmt'].sel(age_cohort = 'age65plus'), _polygons)\n", + "_polygons_mmt = analysis_utils.xarray_to_gpd(\n", + " betas[\"mmt\"].sel(age_cohort=\"age65plus\"), _polygons\n", + ")\n", "\n", "ax = analysis_utils.plot_single(\n", " _polygons_mmt,\n", - " col='mmt',\n", - " cm = 'YlOrRd',\n", - " vmin = 15, vmax = 30, n_colors = 9,\n", - " sup_title=f\"\",\n", - " cbar_label= \"MMT [C]\",\n", + " col=\"mmt\",\n", + " cm=\"YlOrRd\",\n", + " vmin=15,\n", + " vmax=30,\n", + " n_colors=9,\n", + " sup_title=\"\",\n", + " cbar_label=\"MMT [C]\",\n", ")\n", "fig = ax.figure\n", - "fig.savefig(\"2608_mmt.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_mmt.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -1332,18 +1331,20 @@ } ], "source": [ - "_polygons_climtas = analysis_utils.xarray_to_gpd(betas['climtas'], _polygons)\n", + "_polygons_climtas = analysis_utils.xarray_to_gpd(betas[\"climtas\"], _polygons)\n", "\n", "ax = analysis_utils.plot_single(\n", " _polygons_climtas,\n", - " col='climtas',\n", - " cm = 'YlOrRd',\n", - " vmin = 0, vmax = 30, n_colors = 9,\n", - " sup_title=f\"\",\n", - " cbar_label= \"Climtas [C]\",\n", + " col=\"climtas\",\n", + " cm=\"YlOrRd\",\n", + " vmin=0,\n", + " vmax=30,\n", + " n_colors=9,\n", + " sup_title=\"\",\n", + " cbar_label=\"Climtas [C]\",\n", ")\n", "fig = ax.figure\n", - "fig.savefig(\"2608_climtas.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_climtas.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -1372,18 +1373,20 @@ } ], "source": [ - "_polygons_gdppc = analysis_utils.xarray_to_gpd(betas['loggdppc'], _polygons)\n", + "_polygons_gdppc = analysis_utils.xarray_to_gpd(betas[\"loggdppc\"], _polygons)\n", "\n", "ax = analysis_utils.plot_single(\n", " _polygons_gdppc,\n", - " col='loggdppc',\n", - " cm = 'YlOrRd',\n", - " vmin = 5, vmax = 15, n_colors = 10,\n", - " sup_title=f\"\",\n", - " cbar_label= \"GDP Per Cap [2023$]\",\n", + " col=\"loggdppc\",\n", + " cm=\"YlOrRd\",\n", + " vmin=5,\n", + " vmax=15,\n", + " n_colors=10,\n", + " sup_title=\"\",\n", + " cbar_label=\"GDP Per Cap [2023$]\",\n", ")\n", "fig = ax.figure\n", - "fig.savefig(\"2608_loggdppc.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_loggdppc.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -1412,9 +1415,11 @@ "outputs": [], "source": [ "polygons_scatter = _polygons_mmt\n", - "polygons_scatter['mmt'] = polygons_scatter['mmt'].mask(polygons_scatter['mmt'] < 0)\n", - "polygons_scatter['climtas'] = _polygons_climtas['climtas']\n", - "polygons_scatter['gdppc'] = _polygons_gdppc['gdppc'].mask(_polygons_gdppc['loggdppc'] < 0)" + "polygons_scatter[\"mmt\"] = polygons_scatter[\"mmt\"].mask(polygons_scatter[\"mmt\"] < 0)\n", + "polygons_scatter[\"climtas\"] = _polygons_climtas[\"climtas\"]\n", + "polygons_scatter[\"gdppc\"] = _polygons_gdppc[\"gdppc\"].mask(\n", + " _polygons_gdppc[\"loggdppc\"] < 0\n", + ")" ] }, { @@ -1850,7 +1855,7 @@ } ], "source": [ - "polygons_scatter[polygons_scatter['mmt'].isna()]" + "polygons_scatter[polygons_scatter[\"mmt\"].isna()]" ] }, { @@ -2102,7 +2107,7 @@ } ], "source": [ - "polygons_scatter.sort_values('mmt', ascending = True)" + "polygons_scatter.sort_values(\"mmt\", ascending=True)" ] }, { @@ -2123,7 +2128,7 @@ } ], "source": [ - "_polygons_gdppc['loggdppc'].min()" + "_polygons_gdppc[\"loggdppc\"].min()" ] }, { @@ -2154,7 +2159,7 @@ } ], "source": [ - "sns.scatterplot(data=polygons_scatter, x='mmt', y='gdppc', hue='climtas')" + "sns.scatterplot(data=polygons_scatter, x=\"mmt\", y=\"gdppc\", hue=\"climtas\")" ] }, { @@ -2185,7 +2190,7 @@ } ], "source": [ - "sns.scatterplot(data=polygons_scatter, x='mmt', y='climtas', hue = 'gdppc')" + "sns.scatterplot(data=polygons_scatter, x=\"mmt\", y=\"climtas\", hue=\"gdppc\")" ] }, { @@ -2203,7 +2208,9 @@ "metadata": {}, "outputs": [], "source": [ - "gammas = xr.open_zarr('/home/emily_zuetell/projects/poreallas/scripts/data/gammas_local.zarr')" + "gammas = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/scripts/data/gammas_local.zarr\"\n", + ")" ] }, { @@ -2224,7 +2231,7 @@ } ], "source": [ - "gammas['gamma_mean'].sel(age_cohort = 'age65plus', degree = 1, covarname = 'climtas').values" + "gammas[\"gamma_mean\"].sel(age_cohort=\"age65plus\", degree=1, covarname=\"climtas\").values" ] }, { diff --git a/analysis/analysis_utils.py b/analysis/analysis_utils.py index 33cd29a..de3d79f 100644 --- a/analysis/analysis_utils.py +++ b/analysis/analysis_utils.py @@ -1,4 +1,3 @@ - ### Utilities to compute age-weighted mortality impacts relative to a base period # ### Emily Zuetell ### July 7, 2026 @@ -10,49 +9,64 @@ import matplotlib.pyplot as plt import matplotlib.colors as mcolors -def get_baseline_period(effect_xr, years = 30): - #Baseline + +def get_baseline_period(effect_xr, years=30): + # Baseline _max_year = effect_xr["/baseline"]["time"].max().dt.year.item() # Subtracting 29 even though we want 30 year baseline because the time slice is an inclusive range. - _min_year = _max_year - (years-1) + _min_year = _max_year - (years - 1) baseline_period = slice(str(_min_year), str(_max_year)) return baseline_period -def pop_weight_sum(da, socioeconomics, rate=True, age_weight=True, cohort = 'age65plus', impact = True): +def pop_weight_sum( + da, socioeconomics, rate=True, age_weight=True, cohort="age65plus", impact=True +): if age_weight: # Weighted by cohort_population pop_weight = xr.concat( - [socioeconomics['pop0to4'], socioeconomics['pop5to64'], socioeconomics['pop65plus']], - dim=pd.Index(['age0to4', 'age5to64', 'age65plus'], name='age_cohort') + [ + socioeconomics["pop0to4"], + socioeconomics["pop5to64"], + socioeconomics["pop65plus"], + ], + dim=pd.Index(["age0to4", "age5to64", "age65plus"], name="age_cohort"), ) # Total Mortality = rate*pop/100,000 age_weighted_total = da * pop_weight / 100000 if rate: - # Mortality Rate = Total Age-Weighted Mortality/Total Pop (deaths/100k) - age_weighted_total = age_weighted_total * 100000 / socioeconomics['pop'] - #Sum across age cohorts for each Impact Region - regional_sum = age_weighted_total.sum(dim='age_cohort') + # Mortality Rate = Total Age-Weighted Mortality/Total Pop (deaths/100k) + age_weighted_total = age_weighted_total * 100000 / socioeconomics["pop"] + # Sum across age cohorts for each Impact Region + regional_sum = age_weighted_total.sum(dim="age_cohort") regional_sum.name = "age_weighted_impact" if impact else "age_weighted_effect" else: # Return individual cohort regional_sum = da.sel(age_cohort=cohort) if not rate: - # Total Mortality + # Total Mortality _cohortstem = cohort[3:] _col = f"pop{_cohortstem}" - regional_sum = da.sel(age_cohort=cohort)*socioeconomics[_col] + regional_sum = da.sel(age_cohort=cohort) * socioeconomics[_col] regional_sum.name = f"{cohort}_impact" if impact else f"{cohort}_effect" - + return regional_sum -def compute_impact(projected, socioeconomics, baseline_period, - chunks = {'number':-1, 'sample': -1, 'region': 'auto'}, - ensemble = False, hotonly = False, - rate = False, age_weight = True, cohort = 'age65plus'): + +def compute_impact( + projected, + socioeconomics, + baseline_period, + chunks={"number": -1, "sample": -1, "region": "auto"}, + ensemble=False, + hotonly=False, + rate=False, + age_weight=True, + cohort="age65plus", +): """ Calculates the difference between forecast and baseline period mortality effects (deaths/100k or total deaths) (monthly climatology), then applies population-weighted for age cohorts. @@ -85,10 +99,11 @@ def compute_impact(projected, socioeconomics, baseline_period, by pop_weight_sum. """ - if hotonly == 'hotonly': - #Hotonly + if hotonly == "hotonly": + # Hotonly _baseline = ( - projected["/baseline_hotonly"]["effect"].chunk({'region': 'auto'}) + projected["/baseline_hotonly"]["effect"] + .chunk({"region": "auto"}) .sel(time=baseline_period) .groupby("time.month") .mean() @@ -97,18 +112,22 @@ def compute_impact(projected, socioeconomics, baseline_period, _forecast = ( projected["/forecast_hotonly"]["effect"] .chunk(chunks) - .groupby("time.month").mean() + .groupby("time.month") + .mean() ) else: _forecast = ( projected["/forecast_hotonly"]["effect"] .chunk(chunks) - .mean(dim="number").groupby("time.month").mean() + .mean(dim="number") + .groupby("time.month") + .mean() ) - elif hotonly == 'coldonly': - #Coldonly + elif hotonly == "coldonly": + # Coldonly _baseline = ( - projected["/baseline_coldonly"]["effect"].chunk({'region': 'auto'}) + projected["/baseline_coldonly"]["effect"] + .chunk({"region": "auto"}) .sel(time=baseline_period) .groupby("time.month") .mean() @@ -117,35 +136,41 @@ def compute_impact(projected, socioeconomics, baseline_period, _forecast = ( projected["/forecast_coldonly"]["effect"] .chunk(chunks) - .groupby("time.month").mean() + .groupby("time.month") + .mean() ) else: _forecast = ( projected["/forecast_coldonly"]["effect"] .chunk(chunks) - .mean(dim="number").groupby("time.month").mean() + .mean(dim="number") + .groupby("time.month") + .mean() ) else: ## Net _baseline = ( projected["/baseline"]["effect"] - .chunk({'region': 'auto'}) + .chunk({"region": "auto"}) .sel(time=baseline_period) .groupby("time.month") .mean() ) if ensemble: - #Maintain ensmble dimension + # Maintain ensmble dimension _forecast = ( projected["/forecast"]["effect"] .chunk(chunks) - .groupby("time.month").mean() + .groupby("time.month") + .mean() ) else: _forecast = ( projected["/forecast"]["effect"] .chunk(chunks) - .mean(dim="number").groupby("time.month").mean() + .mean(dim="number") + .groupby("time.month") + .mean() ) # Compute Impact @@ -154,25 +179,35 @@ def compute_impact(projected, socioeconomics, baseline_period, impact.name = "impact" impact.attrs["long_name"] = "Temperature mortality impact" impact.attrs["units"] = "Deaths per 100,000 people" if rate else "Deaths" - impact.attrs['hotonly'] = hotonly + impact.attrs["hotonly"] = hotonly # Apply population weighting - regional_sum = pop_weight_sum(impact, socioeconomics, rate=rate, age_weight=age_weight, cohort = cohort, impact = True) + regional_sum = pop_weight_sum( + impact, + socioeconomics, + rate=rate, + age_weight=age_weight, + cohort=cohort, + impact=True, + ) return regional_sum + ### Analysis Functions ### -def xarray_to_gpd(data, polygons, crs = 'ESRI:54030'): + +def xarray_to_gpd(data, polygons, crs="ESRI:54030"): _polygons_data = polygons.merge( data.to_dataframe(name=data.name or "value").reset_index(), on="region", - ) + ) _polygons_data = _polygons_data.cx[:, -60:90].to_crs(crs) return _polygons_data -def compute_stats(da, dim = 'number', polygon = None): + +def compute_stats(da, dim="number", polygon=None): """ From an xarray, return statistics along dimension, 'dim' @@ -180,70 +215,109 @@ def compute_stats(da, dim = 'number', polygon = None): """ mean = da.mean(dim=dim) - std = da.std(dim = dim) - min = da.min(dim = dim) - max = da.max(dim = dim) - p = da.quantile([0.10, 0.17, 0.5, 0.83, 0.9], dim= dim) - - p10 = p.sel(quantile = 0.10, drop = True) - p17 = p.sel(quantile = 0.17, drop = True) - p50 = p.sel(quantile = 0.50, drop = True) - p83 = p.sel(quantile = 0.83, drop = True) - p90 = p.sel(quantile = 0.90, drop = True) + std = da.std(dim=dim) + min = da.min(dim=dim) + max = da.max(dim=dim) + p = da.quantile([0.10, 0.17, 0.5, 0.83, 0.9], dim=dim) + + p10 = p.sel(quantile=0.10, drop=True) + p17 = p.sel(quantile=0.17, drop=True) + p50 = p.sel(quantile=0.50, drop=True) + p83 = p.sel(quantile=0.83, drop=True) + p90 = p.sel(quantile=0.90, drop=True) likely_range = p83 - p17 - ds_out = xr.Dataset({ - "median" : p50, - "p17": p17, - "p83": p83, - "likely_range_IPCC": likely_range, - "mean": mean, - "std": std, - "min": min, - "max": max, - "p10": p10, - "p90": p90, - }) + ds_out = xr.Dataset( + { + "median": p50, + "p17": p17, + "p83": p83, + "likely_range_IPCC": likely_range, + "mean": mean, + "std": std, + "min": min, + "max": max, + "p10": p10, + "p90": p90, + } + ) if polygon is not None: - _polygons_num = polygon.merge(ds_out - .to_dataframe() - .reset_index(), - on="region",) + _polygons_num = polygon.merge( + ds_out.to_dataframe().reset_index(), + on="region", + ) return _polygons_num return ds_out + ###Output Functions ### -def make_csv(effect, socioeconomics, polygon, baseline_period, ensemble = True, hotonly = "net", rate = False, age_weight = True): +def make_csv( + effect, + socioeconomics, + polygon, + baseline_period, + ensemble=True, + hotonly="net", + rate=False, + age_weight=True, +): rate_l = "rate" if rate else "total" - #Compute Impact from Effect - impact = compute_impact(effect, socioeconomics, - ensemble = ensemble, baseline_period=baseline_period, - hotonly = hotonly, rate = rate, - age_weight = age_weight) + # Compute Impact from Effect + impact = compute_impact( + effect, + socioeconomics, + ensemble=ensemble, + baseline_period=baseline_period, + hotonly=hotonly, + rate=rate, + age_weight=age_weight, + ) # Monthly Stats - stat_cols = ['median', 'p17', 'p83', 'likely_range_IPCC', - 'mean','std', - 'min', 'max', - 'p10', 'p90'] - _polygons_impact = compute_stats(impact, dim=["number", "sample"], polygon = polygon) - wide = _polygons_impact.pivot(index=['region', 'ISO'], columns="month", values=stat_cols) + stat_cols = [ + "median", + "p17", + "p83", + "likely_range_IPCC", + "mean", + "std", + "min", + "max", + "p10", + "p90", + ] + _polygons_impact = compute_stats(impact, dim=["number", "sample"], polygon=polygon) + wide = _polygons_impact.pivot( + index=["region", "ISO"], columns="month", values=stat_cols + ) wide.columns = [f"month {m} {stat}" for stat, m in wide.columns] wide = wide.reset_index() - wide.to_csv(f"2608_{hotonly}_{rate_l}_all_stats.csv", index = False) + wide.to_csv(f"2608_{hotonly}_{rate_l}_all_stats.csv", index=False) # 6-month stats - mo6 = impact.sum(dim = 'month') + mo6 = impact.sum(dim="month") _polygons_mo6 = compute_stats(mo6, dim=["number", "sample"], polygon=polygon) - mo6_out = _polygons_mo6[["region", "ISO", - 'median', 'p17', 'p83', 'likely_range_IPCC', - 'mean','std', - 'min', 'max', - 'p10', 'p90']] - mo6_out.to_csv(f"2608_{hotonly}_6mo_{rate_l}_all_stats.csv", index = False) + mo6_out = _polygons_mo6[ + [ + "region", + "ISO", + "median", + "p17", + "p83", + "likely_range_IPCC", + "mean", + "std", + "min", + "max", + "p10", + "p90", + ] + ] + mo6_out.to_csv(f"2608_{hotonly}_6mo_{rate_l}_all_stats.csv", index=False) return + def build_stats_text(da, dim=None, fmt="{:.2f}"): stats = { "mean": float(da.mean(dim=dim)), @@ -271,16 +345,20 @@ def add_stats_annotation(text, ax, loc="upper left"): va=va, fontsize=9, ) + + ##### Plotting Functions ##### from functools import lru_cache import regionmask -#Get and store land data + +# Get and store land data @lru_cache(maxsize=None) def _get_land(crs): - land = gpd.read_file(geodatasets.get_path('naturalearth land')) + land = gpd.read_file(geodatasets.get_path("naturalearth land")) return land.cx[:, -60:90].to_crs(crs) + @lru_cache(maxsize=None) def _get_land_mask(lon_key, lat_key): dummy = xr.DataArray( @@ -290,11 +368,13 @@ def _get_land_mask(lon_key, lat_key): ) return regionmask.defined_regions.natural_earth_v5_0_0.land_110.mask(dummy) + def land_only(da): da = da.rename({"longitude": "lon", "latitude": "lat"}) mask = _get_land_mask(tuple(da.lon.values), tuple(da.lat.values)) return da.where(mask.notnull() & (da.lat > -60)) + def get_step(target_range, min_bins=2, max_bins=8): magnitude = 10 ** np.floor(np.log10(target_range)) for _ in range(3): @@ -306,16 +386,19 @@ def get_step(target_range, min_bins=2, max_bins=8): magnitude /= 10 return magnitude * 10 + def nice_step(target): """Snap an arbitrary step size to the nearest round value (1, 2, 5, or 10 x 10**n).""" magnitude = 10 ** np.floor(np.log10(target)) options = np.array([1, 2, 5, 10]) * magnitude return options[np.argmin(np.abs(options - target))] + def round_bounds(bounds, step): """Remove floating-point drift so bounds land exactly on multiples of step/2.""" return np.round(bounds * 2 / step) * step / 2 + def get_ticks(bounds, step=None, max_ticks=10, symmetric=False): if symmetric: absmax = bounds.max() @@ -327,12 +410,22 @@ def get_ticks(bounds, step=None, max_ticks=10, symmetric=False): stride = int(np.ceil(len(bounds) / max_ticks)) return bounds[::stride] -def make_cmap(bounds, cm='bwr'): + +def make_cmap(bounds, cm="bwr"): cmap_base = plt.get_cmap(cm) colors = cmap_base(np.linspace(0, 1, len(bounds) - 1)) - diverging_cmaps = {'bwr', 'seismic', 'coolwarm', 'RdBu', 'RdYlBu', 'PiYG', 'PRGn', 'BrBG'} - base_name = cm[:-2] if cm.endswith('_r') else cm + diverging_cmaps = { + "bwr", + "seismic", + "coolwarm", + "RdBu", + "RdYlBu", + "PiYG", + "PRGn", + "BrBG", + } + base_name = cm[:-2] if cm.endswith("_r") else cm if base_name in diverging_cmaps: center_idx = len(colors) // 2 colors[center_idx] = [0.95, 0.95, 0.95, 1] @@ -342,10 +435,20 @@ def make_cmap(bounds, cm='bwr'): sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) return cmap, norm, sm -def build_colormap(gdf=None, col=None, cm='bwr', vmin=None, vmax=None, n_colors=None): + +def build_colormap(gdf=None, col=None, cm="bwr", vmin=None, vmax=None, n_colors=None): # Build discrete colormap, branching on diverging vs sequential - diverging_cmaps = {'bwr', 'seismic', 'coolwarm', 'RdBu', 'RdYlBu', 'PiYG', 'PRGn', 'BrBG'} - base_name = cm[:-2] if cm.endswith('_r') else cm + diverging_cmaps = { + "bwr", + "seismic", + "coolwarm", + "RdBu", + "RdYlBu", + "PiYG", + "PRGn", + "BrBG", + } + base_name = cm[:-2] if cm.endswith("_r") else cm if base_name in diverging_cmaps: # Symmetric around zero, with a white center band @@ -355,7 +458,11 @@ def build_colormap(gdf=None, col=None, cm='bwr', vmin=None, vmax=None, n_colors= absmax = math.ceil(gdf[col].abs().quantile(0.95)) step = nice_step((2 * absmax) / n_colors) if n_colors else get_step(2 * absmax) - bounds = np.arange(-np.ceil(absmax / step) * step - step / 2, np.ceil(absmax / step) * step + step, step) + bounds = np.arange( + -np.ceil(absmax / step) * step - step / 2, + np.ceil(absmax / step) * step + step, + step, + ) bounds = round_bounds(bounds, step) cmap, norm, sm = make_cmap(bounds, cm=cm) else: @@ -367,7 +474,9 @@ def build_colormap(gdf=None, col=None, cm='bwr', vmin=None, vmax=None, n_colors= hi = vmax if vmax is not None else gdf[col].max() step = nice_step((hi - lo) / n_colors) if n_colors else get_step(hi - lo) - bounds = np.arange(np.floor(lo / step) * step, np.ceil(hi / step) * step + step, step) + bounds = np.arange( + np.floor(lo / step) * step, np.ceil(hi / step) * step + step, step + ) bounds = round_bounds(bounds, step) cmap = plt.get_cmap(cm, len(bounds) - 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) @@ -376,15 +485,31 @@ def build_colormap(gdf=None, col=None, cm='bwr', vmin=None, vmax=None, n_colors= ticks = get_ticks(bounds, step=step, symmetric=(base_name in diverging_cmaps)) return cmap, norm, sm, ticks, step + import math import geodatasets -land = gpd.read_file(geodatasets.get_path('naturalearth land')) -def plot_single(gdf, col, sup_title="", save_title="", cm='bwr', - cbar_label=None, vmin=None, vmax=None, edgecolor=None, linewidth=0, - ax=None, cbar_location='right', colorbar=True, n_colors=None, - annotation = None, - target_crs = 'ESRI:54030'): +land = gpd.read_file(geodatasets.get_path("naturalearth land")) + + +def plot_single( + gdf, + col, + sup_title="", + save_title="", + cm="bwr", + cbar_label=None, + vmin=None, + vmax=None, + edgecolor=None, + linewidth=0, + ax=None, + cbar_location="right", + colorbar=True, + n_colors=None, + annotation=None, + target_crs="ESRI:54030", +): if gdf.crs is None: raise ValueError("gdf has no CRS set") @@ -392,17 +517,19 @@ def plot_single(gdf, col, sup_title="", save_title="", cm='bwr', gdf = gdf.to_crs(target_crs) land = _get_land(target_crs) - - cmap, norm, sm, ticks, step = build_colormap(gdf, col, cm=cm, vmin=vmin, vmax=vmax, n_colors=n_colors) + + cmap, norm, sm, ticks, step = build_colormap( + gdf, col, cm=cm, vmin=vmin, vmax=vmax, n_colors=n_colors + ) standalone = ax is None if standalone: fig, ax = plt.subplots(figsize=(8, 6)) else: fig = ax.figure - #Plot background land - land.plot(color='lightgray', edgecolor='lightgray', ax=ax) - #Plot data + # Plot background land + land.plot(color="lightgray", edgecolor="lightgray", ax=ax) + # Plot data gdf.plot( column=col, legend=False, @@ -415,73 +542,135 @@ def plot_single(gdf, col, sup_title="", save_title="", cm='bwr', ax.set_axis_off() if colorbar: - orientation = 'horizontal' if cbar_location in ('bottom', 'top') else 'vertical' - cb = fig.colorbar(sm, ax=ax, location=cbar_location, orientation=orientation, shrink=0.6, ticks=ticks, label=cbar_label) + orientation = "horizontal" if cbar_location in ("bottom", "top") else "vertical" + cb = fig.colorbar( + sm, + ax=ax, + location=cbar_location, + orientation=orientation, + shrink=0.6, + ticks=ticks, + label=cbar_label, + ) decimals = max(0, int(np.ceil(-np.log10(step)))) if step < 1 else 0 cb.set_ticklabels([f"{t:.{decimals}f}" for t in ticks]) for label in cb.ax.get_xticklabels(): label.set_rotation(45) - label.set_ha('right') + label.set_ha("right") if annotation: add_stats_annotation(annotation, ax) - + if standalone and save_title: fig.savefig(save_title, dpi=600, bbox_inches="tight") return ax -def plot_monthly(gdf, col, sup_title="", save_title="", - cm='bwr', cbar_label=None, vmin=None, vmax=None, - edgecolor=None, n_colors=None, linewidth=0, - month_order=None): + +def plot_monthly( + gdf, + col, + sup_title="", + save_title="", + cm="bwr", + cbar_label=None, + vmin=None, + vmax=None, + edgecolor=None, + n_colors=None, + linewidth=0, + month_order=None, +): if vmin is None: vmin = -gdf[col].abs().quantile(0.95) if vmax is None: vmax = gdf[col].abs().quantile(0.95) if month_order is None: - months = sorted(gdf['month'].unique()) + months = sorted(gdf["month"].unique()) else: - present = set(gdf['month'].unique()) + present = set(gdf["month"].unique()) months = [m for m in month_order if m in present] fig, axes = plt.subplots(2, 3, figsize=(16, 6)) for ax, month in zip(axes.flat, months): - group = gdf[gdf['month'] == month] + group = gdf[gdf["month"] == month] plot_single( - group, col, cm=cm, vmin=vmin, vmax=vmax, - edgecolor=edgecolor, linewidth=linewidth, ax=ax, colorbar=False, n_colors=n_colors, + group, + col, + cm=cm, + vmin=vmin, + vmax=vmax, + edgecolor=edgecolor, + linewidth=linewidth, + ax=ax, + colorbar=False, + n_colors=n_colors, ) - ax.set_title(f'Month {month}') + ax.set_title(f"Month {month}") - for ax in axes.flat[len(months):]: + for ax in axes.flat[len(months) :]: ax.set_axis_off() - cmap, norm, sm, ticks, step = build_colormap(gdf, col, cm=cm, vmin=vmin, vmax=vmax, n_colors=n_colors) - fig.colorbar(sm, ax=axes.ravel().tolist(), location='right', shrink=0.6, ticks=ticks, label=cbar_label) + cmap, norm, sm, ticks, step = build_colormap( + gdf, col, cm=cm, vmin=vmin, vmax=vmax, n_colors=n_colors + ) + fig.colorbar( + sm, + ax=axes.ravel().tolist(), + location="right", + shrink=0.6, + ticks=ticks, + label=cbar_label, + ) fig.suptitle(sup_title, fontsize=14) fig.savefig(save_title, dpi=600, bbox_inches="tight") return fig -def plot_aggregate(gdf, col, agg='sum', rate=True, sup_title="", save_title="", cm='bwr', - cbar_label=None, vmin=None, vmax=None, edgecolor=None, linewidth=0, - ax=None, cbar_location='right', colorbar=True, n_colors=None, - annotation = None, - target_crs = 'ESRI:54030'): - - annual = gdf.groupby(['region', 'geometry'], as_index=False)[col].agg(agg) - annual = gpd.GeoDataFrame(annual, geometry='geometry', crs=gdf.crs) - - ax = plot_single(annual, col,cm=cm, ax=ax, cbar_label=cbar_label, - vmin=vmin, vmax=vmax, edgecolor=edgecolor, linewidth=linewidth, - cbar_location=cbar_location, colorbar=colorbar,n_colors=n_colors, - annotation=annotation, - target_crs = target_crs) +def plot_aggregate( + gdf, + col, + agg="sum", + rate=True, + sup_title="", + save_title="", + cm="bwr", + cbar_label=None, + vmin=None, + vmax=None, + edgecolor=None, + linewidth=0, + ax=None, + cbar_location="right", + colorbar=True, + n_colors=None, + annotation=None, + target_crs="ESRI:54030", +): + + annual = gdf.groupby(["region", "geometry"], as_index=False)[col].agg(agg) + annual = gpd.GeoDataFrame(annual, geometry="geometry", crs=gdf.crs) + + ax = plot_single( + annual, + col, + cm=cm, + ax=ax, + cbar_label=cbar_label, + vmin=vmin, + vmax=vmax, + edgecolor=edgecolor, + linewidth=linewidth, + cbar_location=cbar_location, + colorbar=colorbar, + n_colors=n_colors, + annotation=annotation, + target_crs=target_crs, + ) fig = ax.figure fig.suptitle(sup_title, fontsize=14) fig.savefig(save_title, dpi=600, bbox_inches="tight") - return fig \ No newline at end of file + return fig diff --git a/analysis/diagnostic_utils.py b/analysis/diagnostic_utils.py index ed3db2b..f0ce8e7 100644 --- a/analysis/diagnostic_utils.py +++ b/analysis/diagnostic_utils.py @@ -1,12 +1,8 @@ import os from dotenv import load_dotenv -import geopandas -import marimo as mo -import matplotlib.pyplot as plt import numpy as np import xarray as xr -from xclim.core import units load_dotenv() @@ -19,56 +15,79 @@ REGIONS_URI = os.environ["POREALLAS_REGIONS_URI"] BETAS_PATH = os.environ["BETAS_PATH"] + def weighted_cdf(data, bins, weights): - + edges = np.concatenate([bins, [bins[-1] + np.diff(bins)[-1]]]) counts, edges = np.histogram(data, bins=edges, density=False) - weighted = counts * weights * (31/counts.sum()) + weighted = counts * weights * (31 / counts.sum()) cdf = np.cumsum(weighted) return bins, cdf -def compute_cumulative_effect(forecast_local, reanalysis_local, region_filter, months, monthly = False, hotonly = False): + +def compute_cumulative_effect( + forecast_local, + reanalysis_local, + region_filter, + months, + monthly=False, + hotonly=False, +): betas_mmt = xr.open_zarr(os.path.join(DATA_DIR, BETAS_PATH)) cdf_data = {} max_cdf = 0 if monthly: for month in months: - forecast_local_month = forecast_local.sel(region = region_filter, time=forecast_local.time.dt.month == month) - reanalysis_local_month = reanalysis_local.sel(region = region_filter, time=reanalysis_local.time.dt.month == month) - - - ref_vals = betas_mmt['mmt'].sel(region=region_filter).sel(age_cohort='age65plus') - betas = betas_mmt['beta_hotonly'] if hotonly else betas_mmt['beta'] - da_temp_bins = betas.sel(region=region_filter).sel(age_cohort='age65plus') - da_temp_bins['tas_bin'] = da_temp_bins['tas_bin'] - - bins = da_temp_bins['tas_bin'].values - weights = betas.sel(region = region_filter).sel(age_cohort = 'age65plus').values - - era5_centers, era5_cdf = weighted_cdf(reanalysis_local_month.values.flatten(), bins, weights) - fc_centers, fc_cdf = weighted_cdf(forecast_local_month.values.flatten(), bins, weights) + forecast_local_month = forecast_local.sel( + region=region_filter, time=forecast_local.time.dt.month == month + ) + reanalysis_local_month = reanalysis_local.sel( + region=region_filter, time=reanalysis_local.time.dt.month == month + ) + + ref_vals = ( + betas_mmt["mmt"].sel(region=region_filter).sel(age_cohort="age65plus") + ) + betas = betas_mmt["beta_hotonly"] if hotonly else betas_mmt["beta"] + da_temp_bins = betas.sel(region=region_filter).sel(age_cohort="age65plus") + da_temp_bins["tas_bin"] = da_temp_bins["tas_bin"] + + bins = da_temp_bins["tas_bin"].values + weights = betas.sel(region=region_filter).sel(age_cohort="age65plus").values + + era5_centers, era5_cdf = weighted_cdf( + reanalysis_local_month.values.flatten(), bins, weights + ) + fc_centers, fc_cdf = weighted_cdf( + forecast_local_month.values.flatten(), bins, weights + ) cdf_data[month] = (era5_centers, era5_cdf, fc_centers, fc_cdf) max_cdf = max(max_cdf, era5_cdf.max(), fc_cdf.max()) else: - forecast_local_month = forecast_local.sel(region = region_filter) - reanalysis_local_month = reanalysis_local.sel(region = region_filter) - - - ref_vals = betas_mmt['mmt'].sel(region=region_filter).sel(age_cohort='age65plus') - betas = betas_mmt['beta_hotonly'] if hotonly else betas_mmt['beta'] - da_temp_bins = betas.sel(region=region_filter).sel(age_cohort='age65plus') - da_temp_bins['tas_bin'] = da_temp_bins['tas_bin'] - - bins = da_temp_bins['tas_bin'].values - weights = betas.sel(region = region_filter).sel(age_cohort = 'age65plus').values - - era5_centers, era5_cdf = weighted_cdf(reanalysis_local_month.values.flatten(), bins, weights) - fc_centers, fc_cdf = weighted_cdf(forecast_local_month.values.flatten(), bins, weights) + forecast_local_month = forecast_local.sel(region=region_filter) + reanalysis_local_month = reanalysis_local.sel(region=region_filter) + + ref_vals = ( + betas_mmt["mmt"].sel(region=region_filter).sel(age_cohort="age65plus") + ) + betas = betas_mmt["beta_hotonly"] if hotonly else betas_mmt["beta"] + da_temp_bins = betas.sel(region=region_filter).sel(age_cohort="age65plus") + da_temp_bins["tas_bin"] = da_temp_bins["tas_bin"] + + bins = da_temp_bins["tas_bin"].values + weights = betas.sel(region=region_filter).sel(age_cohort="age65plus").values + + era5_centers, era5_cdf = weighted_cdf( + reanalysis_local_month.values.flatten(), bins, weights + ) + fc_centers, fc_cdf = weighted_cdf( + forecast_local_month.values.flatten(), bins, weights + ) cdf_data = (era5_centers, era5_cdf, fc_centers, fc_cdf) max_cdf = max(max_cdf, era5_cdf.max(), fc_cdf.max()) - return cdf_data, max_cdf \ No newline at end of file + return cdf_data, max_cdf diff --git a/analysis/isku_utils.py b/analysis/isku_utils.py index 30b7ec8..420d7bf 100644 --- a/analysis/isku_utils.py +++ b/analysis/isku_utils.py @@ -1,13 +1,7 @@ import os from dotenv import load_dotenv -import geopandas -import marimo as mo -import matplotlib.pyplot as plt -import numpy as np -import seaborn as sns import xarray as xr -from xclim.core import units import isku @@ -24,11 +18,13 @@ def _do_nothing(ds: xr.Dataset) -> xr.Dataset: return ds + do_nothing_func = isku.build_extraction_template( pre=_do_nothing, post=lambda ds: ds.astype("float32"), # Save space. Don't need float64. ) + def read_regions(uri: str) -> isku.GridWeightingRegions: _region_weights = xr.load_dataset(uri)[ ["lat", "lon", "region", "weight"] @@ -39,12 +35,13 @@ def read_regions(uri: str) -> isku.GridWeightingRegions: regions = isku.GridWeightingRegions(_region_weights) # ty: ignore[invalid-argument-type] return regions -def grid_to_ir(data, savefile = None): - if 'longitude' in data.dims: - if data['longitude'].min() >= 0: + +def grid_to_ir(data, savefile=None): + if "longitude" in data.dims: + if data["longitude"].min() >= 0: data = lon_adjust(data) else: - data = lon_adjust(data, roll = False) + data = lon_adjust(data, roll=False) regions = read_regions(os.path.join(DATA_DIR, REGIONS_URI)) data_ir = isku.extract_regions( data, @@ -56,7 +53,8 @@ def grid_to_ir(data, savefile = None): data_ir.to_zarr(f"{savefile}.zarr") return data_ir -def lon_adjust(_ds, roll = True): + +def lon_adjust(_ds, roll=True): if roll: _ds["longitude"] = (_ds["longitude"] + 180) % 360 - 180 _ds = _ds.sortby("longitude") diff --git a/analysis/u00_regionalization.ipynb b/analysis/u00_regionalization.ipynb index f122207..0d75772 100644 --- a/analysis/u00_regionalization.ipynb +++ b/analysis/u00_regionalization.ipynb @@ -20,19 +20,17 @@ "source": [ "import os\n", "from pathlib import Path\n", - "from dotenv import load_dotenv, find_dotenv\n", + "from dotenv import load_dotenv\n", "import xarray as xr\n", "\n", "import geopandas as gpd\n", - "import pandas as pd\n", - "import numpy as np\n", - "import math\n", "import matplotlib.pyplot as plt\n", "\n", "import analysis_utils\n", "import isku_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -45,14 +43,10 @@ "outputs": [], "source": [ "import cil_regionalization as cilreg\n", - "from pathlib import Path\n", - "import geopandas as gpd\n", - "import pandas as pd\n", "from cil_regionalization.config import SourceUnitPolicies\n", "\n", "DATA = Path(\"data\")\n", - "DATA_VERSION = \"world-combo-201710\" # the impact region version of the sample\n", - "\n" + "DATA_VERSION = \"world-combo-201710\" # the impact region version of the sample" ] }, { @@ -68,7 +62,7 @@ "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", "\n", - "#Climate Data\n", + "# Climate Data\n", "TAS_FORECAST_URI = os.environ[\"POREALLAS_TAS_FORECAST_URI\"]\n", "ERA5_URI = os.environ[\"POREALLAS_ERA5_URI\"]" ] @@ -80,10 +74,10 @@ "metadata": {}, "outputs": [], "source": [ - "#Projection Effects\n", - "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated = False)\n", - "baseline_period = analysis_utils.get_baseline_period(effect, years = 30)\n", - "#Impact Regions\n", + "# Projection Effects\n", + "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated=False)\n", + "baseline_period = analysis_utils.get_baseline_period(effect, years=30)\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(os.path.join(DATA_DIR, IMPACT_REGION_POLYGONS))\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -91,9 +85,11 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(os.path.join(DATA_DIR, SOCIOECONOMICS_URI))\n", - "socioeconomics = socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop', 'gdppc', 'iso3']]" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(os.path.join(DATA_DIR, SOCIOECONOMICS_URI))\n", + "socioeconomics = socioeconomics.sel(year=2026)[\n", + " [\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\", \"gdppc\", \"iso3\"]\n", + "]" ] }, { @@ -103,23 +99,39 @@ "metadata": {}, "outputs": [], "source": [ - "def ir_to_adm1(effect, polygons, rate = False, hotonly = \"net\", monthly = True, dims = ['number', 'sample'], operation = \"median\"):\n", - " #Compute the Mortality Impact relative to a Baseline Period for the Projected Effects\n", - " impact = analysis_utils.compute_impact(effect, socioeconomics, ensemble = True, baseline_period=baseline_period, hotonly = hotonly, rate = rate, age_weight = True)\n", + "def ir_to_adm1(\n", + " effect,\n", + " polygons,\n", + " rate=False,\n", + " hotonly=\"net\",\n", + " monthly=True,\n", + " dims=[\"number\", \"sample\"],\n", + " operation=\"median\",\n", + "):\n", + " # Compute the Mortality Impact relative to a Baseline Period for the Projected Effects\n", + " impact = analysis_utils.compute_impact(\n", + " effect,\n", + " socioeconomics,\n", + " ensemble=True,\n", + " baseline_period=baseline_period,\n", + " hotonly=hotonly,\n", + " rate=rate,\n", + " age_weight=True,\n", + " )\n", "\n", " if not monthly:\n", - " impact_s = impact.sum(dim = \"month\")\n", + " impact_s = impact.sum(dim=\"month\")\n", " else:\n", " impact_s = impact\n", "\n", " if operation == \"median\":\n", - " impact_s = impact_s.median(dim = dims)\n", + " impact_s = impact_s.median(dim=dims)\n", " elif operation == \"mean\":\n", - " impact_s = impact_s.mean(dim = dims)\n", + " impact_s = impact_s.mean(dim=dims)\n", "\n", " _polygons_impact = analysis_utils.xarray_to_gpd(impact_s, polygons)\n", " shapes = gpd.read_parquet(DATA / \"adm1.parquet\")\n", - " #names = shapes[[\"ISO\", \"ID_1\", \"NAME_1\"]]\n", + " # names = shapes[[\"ISO\", \"ID_1\", \"NAME_1\"]]\n", " names = shapes[[\"GID_0\", \"GID_1\", \"NAME_1\"]]\n", " # Total Deaths: \"per_source\", \"kind = extensive\"\n", " # Death Rates: \"per_destination\", \"kind = intensive\"\n", @@ -127,24 +139,32 @@ " weights = cilreg.fetch_weights(\"gadm41-adm1-per-destination\")\n", " else:\n", " weights = cilreg.fetch_weights(\"gadm41-adm1-per-source\")\n", - " \n", - " damages = _polygons_impact[[\"region\", \"month\", \"age_weighted_impact\"]] if monthly else _polygons_impact[[\"region\", \"age_weighted_impact\"]]\n", - " damages = damages.rename(columns = {\"region\": \"hierid\"})\n", + "\n", + " damages = (\n", + " _polygons_impact[[\"region\", \"month\", \"age_weighted_impact\"]]\n", + " if monthly\n", + " else _polygons_impact[[\"region\", \"age_weighted_impact\"]]\n", + " )\n", + " damages = damages.rename(columns={\"region\": \"hierid\"})\n", " col_regions = {(h,) for h in damages[\"hierid\"].unique()}\n", - " #allow_partial_coverage=True to handle coastal geometry\n", + " # allow_partial_coverage=True to handle coastal geometry\n", " kind = \"intensive\" if rate else \"extensive\"\n", "\n", " _damages = damages\n", " adm1 = cilreg.apply_weights(\n", - " weights, _damages,\n", - " kind=kind, weight=\"pop\", value_col=\"age_weighted_impact\",\n", - " data_version=DATA_VERSION, restrict_to_sources=col_regions,allow_partial_coverage=True,\n", - " policies=SourceUnitPolicies(on_unmatched='skip', on_zero_weight='skip'),\n", + " weights,\n", + " _damages,\n", + " kind=kind,\n", + " weight=\"pop\",\n", + " value_col=\"age_weighted_impact\",\n", + " data_version=DATA_VERSION,\n", + " restrict_to_sources=col_regions,\n", + " allow_partial_coverage=True,\n", + " policies=SourceUnitPolicies(on_unmatched=\"skip\", on_zero_weight=\"skip\"),\n", " ).frame\n", "\n", " m = shapes.merge(adm1, on=[\"GID_0\", \"GID_1\"])\n", "\n", - "\n", " return m" ] }, @@ -163,7 +183,15 @@ } ], "source": [ - "m = ir_to_adm1(effect, _polygons, rate = False, hotonly = \"net\", monthly = False, dims = ['number', 'sample'], operation = \"median\")" + "m = ir_to_adm1(\n", + " effect,\n", + " _polygons,\n", + " rate=False,\n", + " hotonly=\"net\",\n", + " monthly=False,\n", + " dims=[\"number\", \"sample\"],\n", + " operation=\"median\",\n", + ")" ] }, { @@ -173,8 +201,9 @@ "metadata": {}, "outputs": [], "source": [ - "m[['GID_0', 'GID_1', 'NAME_0', 'NAME_1',\n", - " 'age_weighted_impact']].to_csv(\"2608_monthly_hotonly_total_median_adm1.csv\")" + "m[[\"GID_0\", \"GID_1\", \"NAME_0\", \"NAME_1\", \"age_weighted_impact\"]].to_csv(\n", + " \"2608_monthly_hotonly_total_median_adm1.csv\"\n", + ")" ] }, { @@ -205,7 +234,7 @@ } ], "source": [ - "m.plot(column = \"age_weighted_impact\", vmin = -15, vmax = 15, cmap = 'bwr')" + "m.plot(column=\"age_weighted_impact\", vmin=-15, vmax=15, cmap=\"bwr\")" ] }, { @@ -225,9 +254,11 @@ "source": [ "# Compare totals for IR and ADM1\n", "for iso in [\"AUS\"]:\n", - " adm1_iso = m[m[\"GID_0\"]==iso]['age_weighted_impact'].sum()\n", - " ir_iso = _polygons_impact[_polygons_impact['ISO']==iso]['age_weighted_impact'].sum()\n", - " print(f\"{iso}\\nIR Total: {ir_iso}\\n ADM1 Total: {adm1_iso}\")\n" + " adm1_iso = m[m[\"GID_0\"] == iso][\"age_weighted_impact\"].sum()\n", + " ir_iso = _polygons_impact[_polygons_impact[\"ISO\"] == iso][\n", + " \"age_weighted_impact\"\n", + " ].sum()\n", + " print(f\"{iso}\\nIR Total: {ir_iso}\\n ADM1 Total: {adm1_iso}\")" ] }, { @@ -259,10 +290,12 @@ ], "source": [ "# Compare Geometry Outlines\n", - "iso = 'PAK'\n", + "iso = \"PAK\"\n", "fig, ax = plt.subplots()\n", - "_polygons_impact[_polygons_impact['ISO']==iso].to_crs(\"WGS 84\").plot(color = 'blue', alpha = 0.4, edgecolor= 'k', ax = ax)\n", - "m[m['GID_0']==iso].plot(color = 'red', edgecolor = 'k', alpha = 0.4, ax = ax)" + "_polygons_impact[_polygons_impact[\"ISO\"] == iso].to_crs(\"WGS 84\").plot(\n", + " color=\"blue\", alpha=0.4, edgecolor=\"k\", ax=ax\n", + ")\n", + "m[m[\"GID_0\"] == iso].plot(color=\"red\", edgecolor=\"k\", alpha=0.4, ax=ax)" ] }, { @@ -294,28 +327,32 @@ ], "source": [ "# Compare Maps\n", - "fig, axs = plt.subplots(1, 2, figsize = (15, 6))\n", + "fig, axs = plt.subplots(1, 2, figsize=(15, 6))\n", "\n", "analysis_utils.plot_single(\n", " _polygons_impact,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " n_colors = 10,\n", - " vmin = -20, vmax = 20,\n", - " sup_title=f\"\",\n", - " cbar_label= \"deaths/100k\" if rate else \"deaths\",\n", - " cbar_location='bottom', ax = axs[0]\n", + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " n_colors=10,\n", + " vmin=-20,\n", + " vmax=20,\n", + " sup_title=\"\",\n", + " cbar_label=\"deaths/100k\" if rate else \"deaths\",\n", + " cbar_location=\"bottom\",\n", + " ax=axs[0],\n", ")\n", "\n", "analysis_utils.plot_single(\n", " m,\n", - " col='age_weighted_impact',\n", - " cm = 'bwr',\n", - " n_colors = 10,\n", - " vmin = -20, vmax = 20,\n", - " sup_title=f\"\",\n", - " cbar_label= \"deaths/100k\" if rate else \"deaths\",\n", - " cbar_location='bottom', ax = axs[1]\n", + " col=\"age_weighted_impact\",\n", + " cm=\"bwr\",\n", + " n_colors=10,\n", + " vmin=-20,\n", + " vmax=20,\n", + " sup_title=\"\",\n", + " cbar_label=\"deaths/100k\" if rate else \"deaths\",\n", + " cbar_location=\"bottom\",\n", + " ax=axs[1],\n", ")" ] }, diff --git a/analysis/v00_validation_sequence_august.ipynb b/analysis/v00_validation_sequence_august.ipynb index be5a9b3..aad3fcc 100644 --- a/analysis/v00_validation_sequence_august.ipynb +++ b/analysis/v00_validation_sequence_august.ipynb @@ -39,20 +39,17 @@ ], "source": [ "import os\n", - "from pathlib import Path\n", - "from dotenv import load_dotenv, find_dotenv\n", + "from dotenv import load_dotenv\n", "import xarray as xr\n", "\n", "import geopandas as gpd\n", - "import pandas as pd\n", - "import numpy as np\n", - "import math\n", "import matplotlib.pyplot as plt\n", "\n", "import analysis_utils\n", "import isku_utils\n", "\n", "import importlib\n", + "\n", "importlib.reload(analysis_utils)\n", "importlib.reload(isku_utils)" ] @@ -66,12 +63,14 @@ "source": [ "load_dotenv()\n", "DATA_DIR = os.environ[\"DATA_DIR\"]\n", - "#EFFECTS_URI = os.environ[\"POREALLAS_EFFECTS_URI\"]\n", - "EFFECTS_URI = \"/home/emily_zuetell/projects/poreallas/data/2608_effect_gamma_sampled.zarr\"\n", + "# EFFECTS_URI = os.environ[\"POREALLAS_EFFECTS_URI\"]\n", + "EFFECTS_URI = (\n", + " \"/home/emily_zuetell/projects/poreallas/data/2608_effect_gamma_sampled.zarr\"\n", + ")\n", "IMPACT_REGION_POLYGONS = os.environ[\"POREALLAS_REGIONS_POLYGONS_URI\"]\n", "SOCIOECONOMICS_URI = os.environ[\"POREALLAS_SOCIOECONOMICS_URI\"]\n", "\n", - "#Climate Data\n", + "# Climate Data\n", "TAS_FORECAST_URI = os.environ[\"POREALLAS_TAS_FORECAST_URI\"]\n", "ERA5_URI = os.environ[\"POREALLAS_ERA5_URI\"]" ] @@ -83,10 +82,10 @@ "metadata": {}, "outputs": [], "source": [ - "#Projection Effects\n", - "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated = False)\n", - "baseline_period = analysis_utils.get_baseline_period(effect, years = 30)\n", - "#Impact Regions\n", + "# Projection Effects\n", + "effect = xr.open_datatree(os.path.join(DATA_DIR, EFFECTS_URI), consolidated=False)\n", + "baseline_period = analysis_utils.get_baseline_period(effect, years=30)\n", + "# Impact Regions\n", "_polygons = (\n", " gpd.read_parquet(IMPACT_REGION_POLYGONS)\n", " .rename(columns={\"hierid\": \"region\"})\n", @@ -94,9 +93,11 @@ " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "\n", - "#Socioeconomics\n", - "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", - "socioeconomics = socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop', 'gdppc', 'iso3']]" + "# Socioeconomics\n", + "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)\n", + "socioeconomics = socioeconomics.sel(year=2026)[\n", + " [\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\", \"gdppc\", \"iso3\"]\n", + "]" ] }, { @@ -122,10 +123,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Raw Forecast/Hindcast\n", - "s51 = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/parsed/08_ecmwf_parsed.zarr\")\n", - "forecast_raw = s51.sel(time=slice('2026-08-01', '2027-01-31')).groupby(\"time.month\").mean()\n", - "hindcast_raw = s51.sel(time=slice('1993-08-01', '2016-01-31')).groupby(\"time.month\").mean().sel(month = [1, 8, 9, 10, 11, 12])" + "# Raw Forecast/Hindcast\n", + "s51 = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/parsed/08_ecmwf_parsed.zarr\"\n", + ")\n", + "forecast_raw = (\n", + " s51.sel(time=slice(\"2026-08-01\", \"2027-01-31\")).groupby(\"time.month\").mean()\n", + ")\n", + "hindcast_raw = (\n", + " s51.sel(time=slice(\"1993-08-01\", \"2016-01-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + " .sel(month=[1, 8, 9, 10, 11, 12])\n", + ")" ] }, { @@ -135,8 +145,8 @@ "metadata": {}, "outputs": [], "source": [ - "#Compute monthly average anomaly\n", - "anomaly = forecast_raw['tas']-hindcast_raw['tas']\n", + "# Compute monthly average anomaly\n", + "anomaly = forecast_raw[\"tas\"] - hindcast_raw[\"tas\"]\n", "anomaly.name = \"forecast_anomaly\"" ] }, @@ -147,7 +157,7 @@ "metadata": {}, "outputs": [], "source": [ - "anomaly_grid_land = analysis_utils.land_only(anomaly.sel(month = 8).mean(dim = 'number'))" + "anomaly_grid_land = analysis_utils.land_only(anomaly.sel(month=8).mean(dim=\"number\"))" ] }, { @@ -194,10 +204,13 @@ ], "source": [ "## Gridded to Impact Region\n", - "anomaly_ir = isku_utils.grid_to_ir(anomaly.compute(), savefile = None)\n", - "_polygons_anomaly = analysis_utils.xarray_to_gpd(anomaly_ir.mean(dim = 'number'), _polygons)\n", - "_polygons_anomaly_std = analysis_utils.xarray_to_gpd(anomaly_ir.std(dim = 'number'), _polygons)\n", - "\n" + "anomaly_ir = isku_utils.grid_to_ir(anomaly.compute(), savefile=None)\n", + "_polygons_anomaly = analysis_utils.xarray_to_gpd(\n", + " anomaly_ir.mean(dim=\"number\"), _polygons\n", + ")\n", + "_polygons_anomaly_std = analysis_utils.xarray_to_gpd(\n", + " anomaly_ir.std(dim=\"number\"), _polygons\n", + ")" ] }, { @@ -220,15 +233,15 @@ "source": [ "fig = analysis_utils.plot_monthly(\n", " _polygons_anomaly,\n", - " col='value',\n", - " cm = 'bwr',\n", - " n_colors = 9,\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS Anomaly [C]',\n", - " month_order= [8, 9, 10, 11, 12, 1]\n", + " col=\"value\",\n", + " cm=\"bwr\",\n", + " n_colors=9,\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS Anomaly [C]\",\n", + " month_order=[8, 9, 10, 11, 12, 1],\n", ")\n", "\n", - "fig.savefig(\"2608_monthly_anomaly.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_monthly_anomaly.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -251,16 +264,16 @@ "source": [ "fig = analysis_utils.plot_monthly(\n", " _polygons_anomaly_std,\n", - " col='value',\n", - " cm = 'Grays',\n", - " n_colors = 9,\n", - " vmin = 0,\n", - " sup_title=f\"\",\n", - " cbar_label= 'SD TAS Anomaly [C]',\n", - " month_order= [8, 9, 10, 11, 12, 1]\n", + " col=\"value\",\n", + " cm=\"Grays\",\n", + " n_colors=9,\n", + " vmin=0,\n", + " sup_title=\"\",\n", + " cbar_label=\"SD TAS Anomaly [C]\",\n", + " month_order=[8, 9, 10, 11, 12, 1],\n", ")\n", "\n", - "fig.savefig(\"2608_monthly_anomaly_std.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_monthly_anomaly_std.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -281,21 +294,21 @@ } ], "source": [ - "#Validation Annotation\n", - "annotation = analysis_utils.build_stats_text(anomaly_ir.mean(dim = 'month'))\n", - "#6 month average anomaly\n", + "# Validation Annotation\n", + "annotation = analysis_utils.build_stats_text(anomaly_ir.mean(dim=\"month\"))\n", + "# 6 month average anomaly\n", "fig = analysis_utils.plot_aggregate(\n", " _polygons_anomaly,\n", - " col='value',\n", - " agg = 'mean',\n", - " cm = 'bwr',\n", - " n_colors = 7,\n", + " col=\"value\",\n", + " agg=\"mean\",\n", + " cm=\"bwr\",\n", + " n_colors=7,\n", " annotation=annotation,\n", - " sup_title=f\"\",\n", - " cbar_label= 'TAS Anomaly [C]',\n", + " sup_title=\"\",\n", + " cbar_label=\"TAS Anomaly [C]\",\n", ")\n", "\n", - "fig.savefig(\"2608_cum_anomaly_CDS_baseline.png\", dpi = 600, bbox_inches=\"tight\")" + "fig.savefig(\"2608_cum_anomaly_CDS_baseline.png\", dpi=600, bbox_inches=\"tight\")" ] }, { @@ -313,7 +326,7 @@ "metadata": {}, "outputs": [], "source": [ - "forecast = s51.sel(time=slice('2026-08-01', '2026-08-31'))\n", + "forecast = s51.sel(time=slice(\"2026-08-01\", \"2026-08-31\"))\n", "hindcast = s51.sel(time=s51.time.dt.month == 8)" ] }, @@ -1263,7 +1276,7 @@ } ], "source": [ - "da = ptile['tas'].sel(number = 1, time = \"2026-08-02\")\n", + "da = ptile[\"tas\"].sel(number=1, time=\"2026-08-02\")\n", "da_land = analysis_utils.land_only(da)\n", "da_land.plot()" ] @@ -2200,7 +2213,7 @@ } ], "source": [ - "ptile_land['tas'].sel(time = slice('2026-08-02', '2026-08-31'))" + "ptile_land[\"tas\"].sel(time=slice(\"2026-08-02\", \"2026-08-31\"))" ] }, { @@ -2249,9 +2262,9 @@ } ], "source": [ - "da = ptile_land['tas'].sel(time = slice('2026-08-02', '2026-08-31'))\n", + "da = ptile_land[\"tas\"].sel(time=slice(\"2026-08-02\", \"2026-08-31\"))\n", "da = da.assign_coords(time=da.indexes[\"time\"].to_datetimeindex(unsafe=True))\n", - "da.plot.hist(bins = 50, density = True)\n", + "da.plot.hist(bins=50, density=True)\n", "plt.xlabel(\"Forecast Quantile\")\n", "plt.show" ] @@ -2324,7 +2337,7 @@ } ], "source": [ - "qdm_land.af.sel(quantiles = 0, method = 'nearest').plot(vmin = -5, vmax = 5, cmap = 'bwr')" + "qdm_land.af.sel(quantiles=0, method=\"nearest\").plot(vmin=-5, vmax=5, cmap=\"bwr\")" ] }, { @@ -2365,9 +2378,13 @@ "metadata": {}, "outputs": [], "source": [ - "fc_adj = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/2608_forecast_adj_corrected.zarr\")\n", + "fc_adj = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/2608_forecast_adj_corrected.zarr\"\n", + ")\n", "fc_adj = fc_adj.sortby(\"time\")\n", - "era5_adj = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/data/era5_adj_corrected.zarr\")" + "era5_adj = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/data/era5_adj_corrected.zarr\"\n", + ")" ] }, { @@ -2377,7 +2394,9 @@ "metadata": {}, "outputs": [], "source": [ - "era5_raw = xr.open_zarr(\"/home/emily_zuetell/projects/poreallas/scripts/data/era5_daily_1993_2025\")" + "era5_raw = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/scripts/data/era5_daily_1993_2025\"\n", + ")" ] }, { @@ -2418,19 +2437,23 @@ } ], "source": [ - "#lat, lon = 47.1625, 19.5033\n", + "# lat, lon = 47.1625, 19.5033\n", "lat, lon = 28.6139, 77.2090\n", "\n", "fig, ax = plt.subplots()\n", - "da = era5_raw.sel(lat = lat, lon = lon, method = 'nearest')['tas']\n", + "da = era5_raw.sel(lat=lat, lon=lon, method=\"nearest\")[\"tas\"]\n", "da2 = da.assign_coords(year=da.time.dt.year, dayofyear=da.time.dt.dayofyear)\n", "da2 = da2.set_index(time=[\"year\", \"dayofyear\"]).unstack(\"time\")\n", - "da2.plot.line(x=\"dayofyear\", hue=\"year\", color = 'gray', ax = ax, add_legend=False, alpha=0.5)\n", + "da2.plot.line(\n", + " x=\"dayofyear\", hue=\"year\", color=\"gray\", ax=ax, add_legend=False, alpha=0.5\n", + ")\n", "\n", - "da = era5_adj.sel(latitude = lat, longitude = lon, method = 'nearest')['tas']\n", + "da = era5_adj.sel(latitude=lat, longitude=lon, method=\"nearest\")[\"tas\"]\n", "da2 = da.assign_coords(year=da.time.dt.year, dayofyear=da.time.dt.dayofyear)\n", "da2 = da2.set_index(time=[\"year\", \"dayofyear\"]).unstack(\"time\")\n", - "da2.plot.line(x=\"dayofyear\", hue=\"year\", color = 'red',ax = ax, add_legend=False, alpha=0.5)\n", + "da2.plot.line(\n", + " x=\"dayofyear\", hue=\"year\", color=\"red\", ax=ax, add_legend=False, alpha=0.5\n", + ")\n", "\n", "from matplotlib.lines import Line2D\n", "\n", @@ -2441,8 +2464,14 @@ "ax.legend(handles=handles)\n", "\n", "fig, ax = plt.subplots()\n", - "s51.sel(time=slice('2026-08-01', '2027-01-31')).sel(latitude = lat, longitude = lon, method = 'nearest')['tas'].plot.line(x=\"time\", hue=\"number\", color = 'gray', alpha = 0.4, add_legend=False, ax=ax)\n", - "fc_adj.sel(latitude = lat, longitude = lon, method = 'nearest')['tas'].plot.line(x=\"time\", hue=\"number\", color = 'red', alpha = 0.4, add_legend=False, ax=ax)\n", + "s51.sel(time=slice(\"2026-08-01\", \"2027-01-31\")).sel(\n", + " latitude=lat, longitude=lon, method=\"nearest\"\n", + ")[\"tas\"].plot.line(\n", + " x=\"time\", hue=\"number\", color=\"gray\", alpha=0.4, add_legend=False, ax=ax\n", + ")\n", + "fc_adj.sel(latitude=lat, longitude=lon, method=\"nearest\")[\"tas\"].plot.line(\n", + " x=\"time\", hue=\"number\", color=\"red\", alpha=0.4, add_legend=False, ax=ax\n", + ")\n", "\n", "\n", "handles = [\n", diff --git a/scripts/02-bias_adjust_era5.py b/scripts/02-bias_adjust_era5.py index dcd4dd8..ee2448c 100644 --- a/scripts/02-bias_adjust_era5.py +++ b/scripts/02-bias_adjust_era5.py @@ -21,12 +21,26 @@ UID = str(uuid.uuid4()) START_TIME = datetime.datetime.now(datetime.UTC).isoformat() -gmfd = xr.open_dataset(GMFD_URI, engine = 'zarr', chunks = {}, backend_kwargs = {"storage_options": {"token": "anon"}}) +gmfd = xr.open_dataset( + GMFD_URI, + engine="zarr", + chunks={}, + backend_kwargs={"storage_options": {"token": "anon"}}, +) # Fill extreme values gmfd = gmfd.sortby("latitude").chunk({"latitude": -1, "longitude": 30, "time": -1}) -gmfd = gmfd.where(gmfd["tas"] < 1000).interpolate_na(dim="latitude", method="linear").compute() +gmfd = ( + gmfd.where(gmfd["tas"] < 1000) + .interpolate_na(dim="latitude", method="linear") + .compute() +) -era5 = xr.open_dataset(ERA5_URI, engine = 'zarr', chunks = {}, backend_kwargs = {"storage_options": {"token": "anon"}}) +era5 = xr.open_dataset( + ERA5_URI, + engine="zarr", + chunks={}, + backend_kwargs={"storage_options": {"token": "anon"}}, +) ref = gmfd.sel(time=slice(str(HISTREF_START_YEAR), str(HISTREF_STOP_YEAR))) hist = era5.sel(time=slice(str(HISTREF_START_YEAR), str(HISTREF_STOP_YEAR))) diff --git a/scripts/05_project_effects_diagnostic.ipynb b/scripts/05_project_effects_diagnostic.ipynb index 1965d60..fe3ec61 100644 --- a/scripts/05_project_effects_diagnostic.ipynb +++ b/scripts/05_project_effects_diagnostic.ipynb @@ -42,7 +42,7 @@ "REGIONS_URI = os.environ[\"POREALLAS_REGIONS_URI\"]\n", "\n", "# Output\n", - "#EFFECTS_URI = os.getenv(\"POREALLAS_EFFECTS_URI\")\n", + "# EFFECTS_URI = os.getenv(\"POREALLAS_EFFECTS_URI\")\n", "EFFECTS_URI = \"2605_effects_qdm.zarr\"" ] }, @@ -54,7 +54,12 @@ "outputs": [], "source": [ "def read_reanalysis(uri: str) -> xr.Dataset:\n", - " _ds = xr.load_dataset(uri, engine = 'zarr', chunks = {}, backend_kwargs = {\"storage_options\": {\"token\": \"anon\"}},)\n", + " _ds = xr.load_dataset(\n", + " uri,\n", + " engine=\"zarr\",\n", + " chunks={},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )\n", "\n", " # Clean up longitude. The data goes from longitude 0 to 360. It needs to go -180 to 180 in ascending order.\n", " _ds[\"longitude\"] = (_ds[\"longitude\"] + 180) % 360 - 180\n", @@ -66,7 +71,12 @@ "\n", "\n", "def read_forecast_ensemble(uri: str) -> xr.Dataset:\n", - " _ds = xr.load_dataset(uri, engine = 'zarr', chunks = {}, backend_kwargs = {\"storage_options\": {\"token\": \"anon\"}},)\n", + " _ds = xr.load_dataset(\n", + " uri,\n", + " engine=\"zarr\",\n", + " chunks={},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )\n", "\n", " # Clean up longitude. The data goes from longitude 0 to 360. It needs to go -180 to 180 in ascending order.\n", " _ds[\"longitude\"] = (_ds[\"longitude\"] + 180) % 360 - 180\n", @@ -78,7 +88,9 @@ " # Drop months without required number of obs. Forecast ensemble is for a fixed number of days so we expect to usually trim off the last month of the forecast if it is ragged and missing days beyond a threshold.\n", " # EZ update: time dim handling\n", " if \"time\" not in _ds.coords and \"valid_time\" in _ds.coords:\n", - " _ds = _ds.swap_dims({\"forecast_period\": \"valid_time\"}).rename({\"valid_time\": \"time\"})\n", + " _ds = _ds.swap_dims({\"forecast_period\": \"valid_time\"}).rename(\n", + " {\"valid_time\": \"time\"}\n", + " )\n", " _dt_dim = \"time\"\n", " _n_initial = _ds[_dt_dim].size\n", " _number_obs = _ds[_dt_dim].resample(time=\"ME\").count()\n", @@ -109,9 +121,12 @@ "\n", "\n", "def read_regions(uri: str) -> isku.GridWeightingRegions:\n", - " _region_weights = xr.load_dataset(uri, engine = 'zarr', chunks = {}, backend_kwargs = {\"storage_options\": {\"token\": \"anon\"}},)[\n", - " [\"lat\", \"lon\", \"region\", \"weight\"]\n", - " ] # Load only what we need.\n", + " _region_weights = xr.load_dataset(\n", + " uri,\n", + " engine=\"zarr\",\n", + " chunks={},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )[[\"lat\", \"lon\", \"region\", \"weight\"]] # Load only what we need.\n", " # Apparently in this version of xarray the `.load()` method type-hints it'll return a DataArray instead of a Dataset.\n", " # It is a Dataset (I checked). So telling ty to ignore it.\n", " # # TODO: send bug upstream?\n", @@ -120,11 +135,21 @@ "\n", "\n", "def read_gammas(uri: str) -> xr.Dataset:\n", - " return xr.load_dataset(uri, engine = 'zarr', chunks = {}, backend_kwargs = {\"storage_options\": {\"token\": \"anon\"}},)\n", + " return xr.load_dataset(\n", + " uri,\n", + " engine=\"zarr\",\n", + " chunks={},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )\n", "\n", "\n", "def read_socioeconomics(uri: str) -> xr.Dataset:\n", - " return xr.load_dataset(uri, engine = 'zarr', chunks = {}, backend_kwargs = {\"storage_options\": {\"token\": \"anon\"}},)" + " return xr.load_dataset(\n", + " uri,\n", + " engine=\"zarr\",\n", + " chunks={},\n", + " backend_kwargs={\"storage_options\": {\"token\": \"anon\"}},\n", + " )" ] }, { @@ -136,13 +161,13 @@ "source": [ "# Read Data\n", "reanalysis = read_reanalysis(ERA5_URI)\n", - "print('read_forecast_ensemble')\n", + "print(\"read_forecast_ensemble\")\n", "forecast_ensemble = read_forecast_ensemble(TAS_FORECAST_URI)\n", - "print('read_regions')\n", + "print(\"read_regions\")\n", "regions = read_regions(REGIONS_URI)\n", - "print('read_socioeconomics')\n", + "print(\"read_socioeconomics\")\n", "socioeconomics = read_socioeconomics(SOCIOECONOMICS_URI)\n", - "print('read_gammas')\n", + "print(\"read_gammas\")\n", "gammas = read_gammas(GAMMA_URI)" ] }, @@ -180,7 +205,7 @@ " forecast_ensemble,\n", " template=make_tas_monthly_histogram,\n", " regions=regions,\n", - ")\n" + ")" ] }, { @@ -241,8 +266,16 @@ } ], "source": [ - "hist_tas = histogram_forecast_tas[\"histogram_tas\"].sel(region = \"PAK.2.5.27\").sum(\"number\").groupby(\"time.month\").sum(\"time\")\n", - "hist_tas.plot(x=\"month\", y=\"tas_bin\", cmap=\"RdBu\", xlim=(5,11), ylim=(0, 40), add_colorbar=False)" + "hist_tas = (\n", + " histogram_forecast_tas[\"histogram_tas\"]\n", + " .sel(region=\"PAK.2.5.27\")\n", + " .sum(\"number\")\n", + " .groupby(\"time.month\")\n", + " .sum(\"time\")\n", + ")\n", + "hist_tas.plot(\n", + " x=\"month\", y=\"tas_bin\", cmap=\"RdBu\", xlim=(5, 11), ylim=(0, 40), add_colorbar=False\n", + ")" ] }, { @@ -366,8 +399,7 @@ "projected_forecast[\"effect\"].attrs = {\n", " \"units\": \"deaths per 100,000 people\",\n", " \"long_name\": \"Temperature mortality\",\n", - "}\n", - "\n" + "}" ] }, { @@ -410,12 +442,16 @@ ], "source": [ "fig, ax = plt.subplots()\n", - "mean = projected_forecast[\"effect\"].sel(region = 'CHN.11.109.767', time = '2026-07-01', age_cohort = 'age65plus', sample = 1)\n", - "sampled = projected_forecast[\"effect\"].sel(region = 'CHN.11.109.767', time = '2026-07-01', age_cohort = 'age65plus')\n", + "mean = projected_forecast[\"effect\"].sel(\n", + " region=\"CHN.11.109.767\", time=\"2026-07-01\", age_cohort=\"age65plus\", sample=1\n", + ")\n", + "sampled = projected_forecast[\"effect\"].sel(\n", + " region=\"CHN.11.109.767\", time=\"2026-07-01\", age_cohort=\"age65plus\"\n", + ")\n", "\n", - "sns.histplot(sampled.values.flatten(), kde=True, bins=30, ax = ax, stat = 'density')\n", - "sns.histplot(mean.values.flatten(), kde=True, bins=30, ax = ax, stat = 'density', color = 'k')\n", - "ax.set_xlabel('Effect [deaths/100k]')" + "sns.histplot(sampled.values.flatten(), kde=True, bins=30, ax=ax, stat=\"density\")\n", + "sns.histplot(mean.values.flatten(), kde=True, bins=30, ax=ax, stat=\"density\", color=\"k\")\n", + "ax.set_xlabel(\"Effect [deaths/100k]\")" ] }, { diff --git a/scripts/QC03_Uncertainty_Analysis.ipynb b/scripts/QC03_Uncertainty_Analysis.ipynb index c8e8510..6b376ad 100644 --- a/scripts/QC03_Uncertainty_Analysis.ipynb +++ b/scripts/QC03_Uncertainty_Analysis.ipynb @@ -12,10 +12,8 @@ "\n", "import geopandas\n", "import pandas as pd\n", - "import marimo as mo\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", - "import seaborn as sns\n", "import xarray as xr\n", "import math\n", "import matplotlib.colors as mcolors\n", @@ -51,8 +49,8 @@ "metadata": {}, "outputs": [], "source": [ - "projected = xr.open_datatree(EFFECTS_URI, engine = 'zarr')\n", - "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)" + "projected = xr.open_datatree(EFFECTS_URI, engine=\"zarr\")\n", + "socioeconomics = xr.open_zarr(SOCIOECONOMICS_URI)" ] }, { @@ -63,21 +61,24 @@ "outputs": [], "source": [ "_polygons = (\n", - " geopandas.read_parquet(\"/home/emily_zuetell/projects/poreallas/scripts/data/impact_region.parquet\")\n", + " geopandas.read_parquet(\n", + " \"/home/emily_zuetell/projects/poreallas/scripts/data/impact_region.parquet\"\n", + " )\n", " .rename(columns={\"hierid\": \"region\"})\n", " .set_index(\"region\")\n", " .set_crs(epsg=4326) # Assuming the data is WGS-82.\n", ")\n", "_polygons_effect = _polygons.merge(\n", " projected[\"/forecast\"][\"effect\"]\n", - " .mean(dim=\"number\").groupby(\"time.month\")\n", - " .mean()\n", + " .mean(dim=\"number\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", " .to_dataframe()\n", " .reset_index(),\n", " on=\"region\",\n", ")\n", "_polygons_socioecon = _polygons.merge(\n", - " socioeconomics.sel(year = 2026)[['pop0to4', 'pop5to64', 'pop65plus', 'pop']]\n", + " socioeconomics.sel(year=2026)[[\"pop0to4\", \"pop5to64\", \"pop65plus\", \"pop\"]]\n", " .to_dataframe()\n", " .reset_index(),\n", " on=\"region\",\n", @@ -96,8 +97,8 @@ "metadata": {}, "outputs": [], "source": [ - "da = projected['/forecast_hotonly']['effect']\n", - "da_baseline = projected['/baseline_hotonly']['effect']" + "da = projected[\"/forecast_hotonly\"][\"effect\"]\n", + "da_baseline = projected[\"/baseline_hotonly\"][\"effect\"]" ] }, { @@ -118,32 +119,44 @@ } ], "source": [ - "region_filters = ['THA.76.909.5787', 'PAK.2.5.27', \"BRA.7.804.1862\", \"EGY.8\"] \n", + "region_filters = [\"THA.76.909.5787\", \"PAK.2.5.27\", \"BRA.7.804.1862\", \"EGY.8\"]\n", "\n", - "fig, axes = plt.subplots(1, len(region_filters), figsize=(6*len(region_filters), 5), sharey=True)\n", + "fig, axes = plt.subplots(\n", + " 1, len(region_filters), figsize=(6 * len(region_filters), 5), sharey=True\n", + ")\n", "\n", "for ax, region_filter in zip(axes, region_filters):\n", - " ds = da.sel(region=region_filter, age_cohort='age65plus')\n", - " ds_baseline = da_baseline.sel(region=region_filter, age_cohort='age65plus')\n", + " ds = da.sel(region=region_filter, age_cohort=\"age65plus\")\n", + " ds_baseline = da_baseline.sel(region=region_filter, age_cohort=\"age65plus\")\n", "\n", - " ds = ds.assign_coords(month=ds['time'].dt.month)\n", - " months_present = np.unique(ds['month'].values)\n", + " ds = ds.assign_coords(month=ds[\"time\"].dt.month)\n", + " months_present = np.unique(ds[\"month\"].values)\n", "\n", - " for n in ds['number'].values:\n", - " ds.sel(number=n).plot.line(x='month', add_legend=False, alpha=0.6, linewidth=0.8, color='pink', ax=ax)\n", + " for n in ds[\"number\"].values:\n", + " ds.sel(number=n).plot.line(\n", + " x=\"month\", add_legend=False, alpha=0.6, linewidth=0.8, color=\"pink\", ax=ax\n", + " )\n", "\n", - " for year, da_year in ds_baseline.groupby('time.year'):\n", - " da_year = da_year.assign_coords(month=da_year['time'].dt.month)\n", - " da_year = da_year.sel(time=da_year['month'].isin(months_present))\n", - " da_year.plot.line(x='month', add_legend=False, alpha=0.6, linewidth=0.8, color='gray', ax=ax)\n", + " for year, da_year in ds_baseline.groupby(\"time.year\"):\n", + " da_year = da_year.assign_coords(month=da_year[\"time\"].dt.month)\n", + " da_year = da_year.sel(time=da_year[\"month\"].isin(months_present))\n", + " da_year.plot.line(\n", + " x=\"month\", add_legend=False, alpha=0.6, linewidth=0.8, color=\"gray\", ax=ax\n", + " )\n", "\n", - " ds.mean(dim='number').plot.line(x='month', color='red', linewidth=2, ax=ax, label='Forecast')\n", - " ds_baseline_month = ds_baseline.assign_coords(month=ds_baseline['time'].dt.month)\n", - " ds_baseline_month = ds_baseline_month.sel(time=ds_baseline_month['month'].isin(months_present))\n", - " ds_baseline_month.groupby('month').mean().plot.line(x='month', color='black', linewidth=2, ax=ax, label='Baseline (1995-2025)')\n", + " ds.mean(dim=\"number\").plot.line(\n", + " x=\"month\", color=\"red\", linewidth=2, ax=ax, label=\"Forecast\"\n", + " )\n", + " ds_baseline_month = ds_baseline.assign_coords(month=ds_baseline[\"time\"].dt.month)\n", + " ds_baseline_month = ds_baseline_month.sel(\n", + " time=ds_baseline_month[\"month\"].isin(months_present)\n", + " )\n", + " ds_baseline_month.groupby(\"month\").mean().plot.line(\n", + " x=\"month\", color=\"black\", linewidth=2, ax=ax, label=\"Baseline (1995-2025)\"\n", + " )\n", "\n", " ax.set_xticks(months_present)\n", - " ax.set_xlabel('Month')\n", + " ax.set_xlabel(\"Month\")\n", " ax.set_title(region_filter)\n", " ax.legend()" ] @@ -176,39 +189,51 @@ } ], "source": [ - "region_filters = ['THA.76.909.5787', 'PAK.2.5.27', \"BRA.7.804.1862\", \"EGY.8\"] \n", - "#region_filters = ['THA.76.909.5787', 'PAK.2.5.27'] \n", + "region_filters = [\"THA.76.909.5787\", \"PAK.2.5.27\", \"BRA.7.804.1862\", \"EGY.8\"]\n", + "# region_filters = ['THA.76.909.5787', 'PAK.2.5.27']\n", "\n", - "fig, axes = plt.subplots(1, len(region_filters), figsize=(6*len(region_filters), 5), sharey=True)\n", + "fig, axes = plt.subplots(\n", + " 1, len(region_filters), figsize=(6 * len(region_filters), 5), sharey=True\n", + ")\n", "\n", "offset = 0.15\n", "\n", "for ax, region_filter in zip(axes, region_filters):\n", - " ds = da.sel(region=region_filter, age_cohort='age65plus')\n", - " ds_baseline = da_baseline.sel(region=region_filter, age_cohort='age65plus')\n", + " ds = da.sel(region=region_filter, age_cohort=\"age65plus\")\n", + " ds_baseline = da_baseline.sel(region=region_filter, age_cohort=\"age65plus\")\n", "\n", - " ds = ds.assign_coords(month=ds['time'].dt.month)\n", - " months_present = np.unique(ds['month'].values)\n", + " ds = ds.assign_coords(month=ds[\"time\"].dt.month)\n", + " months_present = np.unique(ds[\"month\"].values)\n", "\n", - " ds_baseline = ds_baseline.assign_coords(month=ds_baseline['time'].dt.month)\n", - " ds_baseline = ds_baseline.sel(time=ds_baseline['month'].isin(months_present))\n", + " ds_baseline = ds_baseline.assign_coords(month=ds_baseline[\"time\"].dt.month)\n", + " ds_baseline = ds_baseline.sel(time=ds_baseline[\"month\"].isin(months_present))\n", "\n", - " data_new = [ds.sel(time=ds['month'] == m).values.flatten() for m in months_present]\n", - " data_baseline = [ds_baseline.sel(time=ds_baseline['month'] == m).values.flatten() for m in months_present]\n", + " data_new = [ds.sel(time=ds[\"month\"] == m).values.flatten() for m in months_present]\n", + " data_baseline = [\n", + " ds_baseline.sel(time=ds_baseline[\"month\"] == m).values.flatten()\n", + " for m in months_present\n", + " ]\n", "\n", - " vp_new = ax.violinplot(data_new, positions=months_present - offset, showmeans=True, widths=0.25)\n", - " vp_baseline = ax.violinplot(data_baseline, positions=months_present + offset, showmeans=True, widths=0.25)\n", + " vp_new = ax.violinplot(\n", + " data_new, positions=months_present - offset, showmeans=True, widths=0.25\n", + " )\n", + " vp_baseline = ax.violinplot(\n", + " data_baseline, positions=months_present + offset, showmeans=True, widths=0.25\n", + " )\n", "\n", - " for body in vp_new['bodies']:\n", - " body.set_facecolor('pink')\n", - " for body in vp_baseline['bodies']:\n", - " body.set_facecolor('gray')\n", + " for body in vp_new[\"bodies\"]:\n", + " body.set_facecolor(\"pink\")\n", + " for body in vp_baseline[\"bodies\"]:\n", + " body.set_facecolor(\"gray\")\n", "\n", " ax.set_xticks(months_present)\n", - " ax.set_xlabel('Month')\n", + " ax.set_xlabel(\"Month\")\n", " ax.set_title(region_filter)\n", "\n", - "axes[0].legend([vp_new['bodies'][0], vp_baseline['bodies'][0]], ['Forecast', 'Baseline (1995-2025)'])" + "axes[0].legend(\n", + " [vp_new[\"bodies\"][0], vp_baseline[\"bodies\"][0]],\n", + " [\"Forecast\", \"Baseline (1995-2025)\"],\n", + ")" ] }, { @@ -218,21 +243,22 @@ "metadata": {}, "outputs": [], "source": [ - "da = projected['/forecast']['effect']\n", + "da = projected[\"/forecast\"][\"effect\"]\n", "p10 = da.quantile(0.10, dim=\"number\")\n", "p90 = da.quantile(0.90, dim=\"number\")\n", "prange = p90 - p10\n", "\n", - "ds_out = xr.Dataset({\n", - " \"p10\": p10.drop_vars(\"quantile\"),\n", - " \"p90\": p90.drop_vars(\"quantile\"),\n", - " \"range\": prange.drop_vars(\"quantile\", errors=\"ignore\")\n", - "})\n", - "_polygons_num = _polygons.merge(ds_out.groupby(\"time.month\")\n", - " .mean()\n", - " .to_dataframe()\n", - " .reset_index(),\n", - " on=\"region\",)" + "ds_out = xr.Dataset(\n", + " {\n", + " \"p10\": p10.drop_vars(\"quantile\"),\n", + " \"p90\": p90.drop_vars(\"quantile\"),\n", + " \"range\": prange.drop_vars(\"quantile\", errors=\"ignore\"),\n", + " }\n", + ")\n", + "_polygons_num = _polygons.merge(\n", + " ds_out.groupby(\"time.month\").mean().to_dataframe().reset_index(),\n", + " on=\"region\",\n", + ")" ] }, { @@ -263,26 +289,27 @@ } ], "source": [ - "months = sorted(_polygons_num['month'].unique())\n", + "months = sorted(_polygons_num[\"month\"].unique())\n", "n = len(months)\n", "ncols = math.ceil(n / 2)\n", - "fig, axes = plt.subplots(4, ncols, figsize=(5*ncols, 14),\n", - " gridspec_kw={'height_ratios': [3, 1, 3, 1]})\n", + "fig, axes = plt.subplots(\n", + " 4, ncols, figsize=(5 * ncols, 14), gridspec_kw={\"height_ratios\": [3, 1, 3, 1]}\n", + ")\n", "\n", "map_axes = list(axes[0]) + list(axes[2])\n", "hist_axes = list(axes[1]) + list(axes[3])\n", "\n", "norm = plt.Normalize(vmin=0, vmax=100)\n", - "sm = plt.cm.ScalarMappable(cmap='viridis', norm=norm)\n", + "sm = plt.cm.ScalarMappable(cmap=\"viridis\", norm=norm)\n", "\n", "xmin, xmax = 0, 10\n", "\n", "for map_ax, hist_ax, m in zip(map_axes, hist_axes, months):\n", - " subset = _polygons_num[_polygons_num['month']==m]\n", - " subset.plot(column='range', ax=map_ax, norm=norm, cmap='viridis', legend=False)\n", + " subset = _polygons_num[_polygons_num[\"month\"] == m]\n", + " subset.plot(column=\"range\", ax=map_ax, norm=norm, cmap=\"viridis\", legend=False)\n", " map_ax.set_title(f\"Month {m}\")\n", "\n", - " hist_ax.hist(subset['range'], bins=100, range=(xmin, xmax))\n", + " hist_ax.hist(subset[\"range\"], bins=100, range=(xmin, xmax))\n", " hist_ax.set_xlim(xmin, xmax)\n", " hist_ax.set_ylim(0, 7000)\n", "\n", @@ -298,8 +325,10 @@ " hist_pos = hist_ax.get_position()\n", " hist_ax.set_position([map_pos.x0, hist_pos.y0, map_pos.width, hist_pos.height])\n", "\n", - "cax = fig.add_axes([0.92, 0.15, 0.015, 0.7]) # dedicated colorbar axis, adjust as needed\n", - "fig.colorbar(sm, cax=cax, label='range')" + "cax = fig.add_axes(\n", + " [0.92, 0.15, 0.015, 0.7]\n", + ") # dedicated colorbar axis, adjust as needed\n", + "fig.colorbar(sm, cax=cax, label=\"range\")" ] }, { @@ -322,14 +351,14 @@ "source": [ "fig, ax = plt.subplots(figsize=(10, 6))\n", "\n", - "months = sorted(_polygons_num['month'].unique())\n", - "data = [_polygons_num.loc[_polygons_num['month'] == m, 'range'] for m in months]\n", + "months = sorted(_polygons_num[\"month\"].unique())\n", + "data = [_polygons_num.loc[_polygons_num[\"month\"] == m, \"range\"] for m in months]\n", "\n", "ax.violinplot(data, positions=months, showmeans=True)\n", "ax.set_xticks(months)\n", - "ax.set_xlabel('Month')\n", - "ax.set_ylabel('10th-90th ptile Range of Mortality Effects')\n", - "ax.set_yscale('log')" + "ax.set_xlabel(\"Month\")\n", + "ax.set_ylabel(\"10th-90th ptile Range of Mortality Effects\")\n", + "ax.set_yscale(\"log\")" ] }, { @@ -347,7 +376,7 @@ "metadata": {}, "outputs": [], "source": [ - "#Hotonly\n", + "# Hotonly\n", "_baseline = (\n", " projected[\"/baseline_hotonly\"][\"effect\"]\n", " .sel(time=baseline_period)\n", @@ -355,13 +384,16 @@ " .mean()\n", ")\n", "_forecast = (\n", - " projected[\"/forecast_hotonly\"][\"effect\"].mean(dim=\"number\").groupby(\"time.month\").mean()\n", + " projected[\"/forecast_hotonly\"][\"effect\"]\n", + " .mean(dim=\"number\")\n", + " .groupby(\"time.month\")\n", + " .mean()\n", ")\n", "impact = _forecast - _baseline\n", "impact.name = \"impact\"\n", "impact.attrs[\"long_name\"] = \"Temperature mortality impact\"\n", "impact.attrs[\"units\"] = \"Deaths per 100,000 people\"\n", - "impact.attrs['hotonly'] = True\n", + "impact.attrs[\"hotonly\"] = True\n", "\n", "_polygons_impact_hotonly = _polygons.merge(\n", " impact.to_dataframe().reset_index(),\n", @@ -381,7 +413,7 @@ "impact.name = \"impact\"\n", "impact.attrs[\"long_name\"] = \"Temperature mortality impact\"\n", "impact.attrs[\"units\"] = \"Deaths per 100,000 people\"\n", - "impact.attrs['hotonly'] = False\n", + "impact.attrs[\"hotonly\"] = False\n", "\n", "_polygons_impact_net = _polygons.merge(\n", " impact.to_dataframe().reset_index(),\n", @@ -405,22 +437,38 @@ "outputs": [], "source": [ "ds_forecast = xr.open_zarr(\"./data/tas_forecast.zarr\")\n", - "ds_forecast = ds_forecast.swap_dims({'forecast_period': 'valid_time'})\n", + "ds_forecast = ds_forecast.swap_dims({\"forecast_period\": \"valid_time\"})\n", "\n", "## Hindcast for Bias Adjustment ##\n", - "ds_hindcast = xr.open_dataset('./data/s51_hcm_allmonths.nc')\n", - "da_hindcast = ds_hindcast['t2m'].assign_coords(forecast_reference_time=ds_hindcast.forecast_reference_time.dt.month).rename({'forecast_reference_time':'month'})\n", + "ds_hindcast = xr.open_dataset(\"./data/s51_hcm_allmonths.nc\")\n", + "da_hindcast = (\n", + " ds_hindcast[\"t2m\"]\n", + " .assign_coords(forecast_reference_time=ds_hindcast.forecast_reference_time.dt.month)\n", + " .rename({\"forecast_reference_time\": \"month\"})\n", + ")\n", "\n", "ds_era5_hcperiod = xr.open_dataset(\"./data/era5_monthly_1981_2016.nc\")\n", - "da_era5 = ds_era5_hcperiod['t2m'].sel(valid_time = slice(\"1993-01-01\", \"2016-12-01\")).groupby(\"valid_time.month\").mean()\n", + "da_era5 = (\n", + " ds_era5_hcperiod[\"t2m\"]\n", + " .sel(valid_time=slice(\"1993-01-01\", \"2016-12-01\"))\n", + " .groupby(\"valid_time.month\")\n", + " .mean()\n", + ")\n", "\n", - "#regridder = xe.Regridder(da_era5, da_hindcast, \"bilinear\")\n", - "#da_era5_regridded = regridder(da_era5)\n", + "# regridder = xe.Regridder(da_era5, da_hindcast, \"bilinear\")\n", + "# da_era5_regridded = regridder(da_era5)\n", "da_era5_regridded = da_era5.interp_like(da_hindcast, method=\"linear\")\n", "\n", "## Historical Baselin ##\n", - "era5_mortality_baseline = xr.open_zarr('/home/emily_zuetell/projects/poreallas/scripts/data/tas_ERA5.zarr')\n", - "era5_mortality_baseline = era5_mortality_baseline['tas'].sel(time = slice(\"1995-01-01\", \"2025-12-31\")).groupby(\"time.month\").mean()" + "era5_mortality_baseline = xr.open_zarr(\n", + " \"/home/emily_zuetell/projects/poreallas/scripts/data/tas_ERA5.zarr\"\n", + ")\n", + "era5_mortality_baseline = (\n", + " era5_mortality_baseline[\"tas\"]\n", + " .sel(time=slice(\"1995-01-01\", \"2025-12-31\"))\n", + " .groupby(\"time.month\")\n", + " .mean()\n", + ")" ] }, { @@ -438,7 +486,6 @@ "metadata": {}, "outputs": [], "source": [ - "import cartopy.crs as ccrs\n", "from scipy.interpolate import interp1d\n", "\n", "rate = False\n", @@ -446,21 +493,23 @@ "age_weight = True\n", "use_log = False\n", "density_hist = True\n", - "regional_sum = pop_weight_sum(col='impact', rate=rate, hotonly = hotonly, age_weight=age_weight)\n", + "regional_sum = pop_weight_sum(\n", + " col=\"impact\", rate=rate, hotonly=hotonly, age_weight=age_weight\n", + ")\n", "\n", - "#region_filter = \"PAK.2.5.27\"\n", - "#region_filter = \"IND.16.203.735\"\n", - "#region_filter = \"BRA.7.804.1862\"\n", - "#region_filter = \"IND.29.439.1743\"\n", + "# region_filter = \"PAK.2.5.27\"\n", + "# region_filter = \"IND.16.203.735\"\n", + "# region_filter = \"BRA.7.804.1862\"\n", + "# region_filter = \"IND.29.439.1743\"\n", "region_filter = \"CHN.11.109.767\" # set to None for global\n", - "#region_filter = \"COD.3.15.51\" \n", - "#region_filter = \"THA.76.909.5787\"\n", - "#region_filter = \"HRV.R4acd3465906d6fba\"\n", - "#region_filter = \"EGY.8\"\n", - "#region_filter = 'BRA.9.R7ae8d9b7a4e0454f'\n", - "#region_filter = 'ZMB.9.67'\n", - "\n", - "months = sorted(regional_sum['month'].unique())[1:]\n", + "# region_filter = \"COD.3.15.51\"\n", + "# region_filter = \"THA.76.909.5787\"\n", + "# region_filter = \"HRV.R4acd3465906d6fba\"\n", + "# region_filter = \"EGY.8\"\n", + "# region_filter = 'BRA.9.R7ae8d9b7a4e0454f'\n", + "# region_filter = 'ZMB.9.67'\n", + "\n", + "months = sorted(regional_sum[\"month\"].unique())[1:]\n", "base_month = months[0]\n", "lead_by_month = {month: (month - base_month + 1) for month in months}\n", "\n", @@ -489,9 +538,14 @@ "hist_axes = []\n", "\n", "for pair_row in range(n_pairs):\n", - " inner_gs = GridSpecFromSubplotSpec(3, ncols, subplot_spec=gs[pair_row],\n", - " hspace=0.15, wspace=0.25,\n", - " height_ratios=[2, 2, 1.2])\n", + " inner_gs = GridSpecFromSubplotSpec(\n", + " 3,\n", + " ncols,\n", + " subplot_spec=gs[pair_row],\n", + " hspace=0.15,\n", + " wspace=0.25,\n", + " height_ratios=[2, 2, 1.2],\n", + " )\n", " for col in range(ncols):\n", " i = pair_row * ncols + col\n", " if i >= n_months:\n", @@ -508,13 +562,12 @@ "bias_adj_iso_all = []\n", "for month in months:\n", " lead = lead_by_month[month]\n", - " bias_adjustment = (\n", - " da_hindcast.sel(month=(month - lead + 1), forecastMonth=lead)\n", - " - da_era5_regridded.sel(month=month)\n", - " )\n", + " bias_adjustment = da_hindcast.sel(\n", + " month=(month - lead + 1), forecastMonth=lead\n", + " ) - da_era5_regridded.sel(month=month)\n", " bias_adj_fc = (\n", - " ds_forecast.sel(valid_time=ds_forecast.valid_time.dt.month == month)['tas']\n", - " .mean(dim=['number', 'valid_time'])\n", + " ds_forecast.sel(valid_time=ds_forecast.valid_time.dt.month == month)[\"tas\"]\n", + " .mean(dim=[\"number\", \"valid_time\"])\n", " .squeeze()\n", " - era5_mortality_baseline.sel(month=month)\n", " - bias_adjustment\n", @@ -522,14 +575,18 @@ " bias_adj_all.append(bias_adj_fc)\n", "\n", " if region_filter:\n", + "\n", " def roll_lon(da):\n", " return da.assign_coords(\n", " longitude=((da.longitude + 180) % 360 - 180)\n", - " ).sortby('longitude')\n", + " ).sortby(\"longitude\")\n", "\n", " bias_adj_fc_iso = (\n", - " roll_lon(ds_forecast.sel(valid_time=ds_forecast.valid_time.dt.month == month)['tas'])\n", - " .squeeze()\n", + " roll_lon(\n", + " ds_forecast.sel(valid_time=ds_forecast.valid_time.dt.month == month)[\n", + " \"tas\"\n", + " ]\n", + " ).squeeze()\n", " - roll_lon(bias_adjustment)\n", " ).sel(\n", " latitude=slice(maxy + bbox_pad, miny - bbox_pad),\n", @@ -540,42 +597,55 @@ "\n", "# --- Shared impact norm ---\n", "\n", - "_polygons_all = _polygons.merge(regional_sum.reset_index(), on=\"region\").to_crs(\"ESRI:54030\")\n", - "absmax_impact = _polygons_all['age_weighted_total'].abs().quantile(0.95)\n", + "_polygons_all = _polygons.merge(regional_sum.reset_index(), on=\"region\").to_crs(\n", + " \"ESRI:54030\"\n", + ")\n", + "absmax_impact = _polygons_all[\"age_weighted_total\"].abs().quantile(0.95)\n", "if use_log:\n", - " norm_impact = mcolors.SymLogNorm(linthresh=1, vmin=-absmax_impact, vmax=absmax_impact)\n", + " norm_impact = mcolors.SymLogNorm(\n", + " linthresh=1, vmin=-absmax_impact, vmax=absmax_impact\n", + " )\n", "else:\n", - " norm_impact = mcolors.TwoSlopeNorm(vmin=-absmax_impact, vcenter=0, vmax=absmax_impact)\n", + " norm_impact = mcolors.TwoSlopeNorm(\n", + " vmin=-absmax_impact, vcenter=0, vmax=absmax_impact\n", + " )\n", "\n", "era5_rolled = era5_daily.assign_coords(\n", " longitude=((era5_daily.longitude + 180) % 360 - 180)\n", - ").sortby('longitude')\n", + ").sortby(\"longitude\")\n", "\n", "era5_local = era5_rolled.sel(\n", " latitude=slice(maxy + bbox_pad, miny - bbox_pad),\n", " longitude=slice(minx - bbox_pad, maxx + bbox_pad),\n", ")\n", "era5_local = (\n", - " era5_local['tas']\n", + " era5_local[\"tas\"]\n", " .sel(time=slice(\"1995-01-01\", \"2025-12-31\"))\n", " .isel(time=era5_local.time.dt.month.isin(months))\n", ")\n", - "xmin = min(era5_local.values.flatten().min(), min(b.values.flatten().min() for b in bias_adj_iso_all))\n", - "xmax = max(era5_local.values.flatten().max(), max(b.values.flatten().max() for b in bias_adj_iso_all))\n", + "xmin = min(\n", + " era5_local.values.flatten().min(),\n", + " min(b.values.flatten().min() for b in bias_adj_iso_all),\n", + ")\n", + "xmax = max(\n", + " era5_local.values.flatten().max(),\n", + " max(b.values.flatten().max() for b in bias_adj_iso_all),\n", + ")\n", "\n", "# --- Plot ---\n", "for i, month in enumerate(months):\n", " ax_top = top_axes[i]\n", " ax_bot = bot_axes[i]\n", " ax_hist = hist_axes[i]\n", - " month_label = pd.Timestamp(2026, month, 1).strftime('%B')\n", + " month_label = pd.Timestamp(2026, month, 1).strftime(\"%B\")\n", "\n", " # Top: temperature\n", " im_temp = bias_adj_all[i].plot(\n", " ax=ax_top,\n", " transform=ccrs.PlateCarree(),\n", - " vmin=-5, vmax=5,\n", - " cmap='RdBu_r',\n", + " vmin=-5,\n", + " vmax=5,\n", + " cmap=\"RdBu_r\",\n", " add_colorbar=False,\n", " )\n", " ax_top.set_title(month_label)\n", @@ -583,12 +653,12 @@ " if region_filter:\n", " ax_top.set_extent(\n", " [minx - bbox_pad, maxx + bbox_pad, miny - bbox_pad, maxy + bbox_pad],\n", - " crs=ccrs.PlateCarree()\n", + " crs=ccrs.PlateCarree(),\n", " )\n", "\n", " # Middle: impact\n", " _polygons_data = _polygons.merge(\n", - " regional_sum[regional_sum['month'] == month].reset_index(),\n", + " regional_sum[regional_sum[\"month\"] == month].reset_index(),\n", " on=\"region\",\n", " )\n", " _plot_data = _polygons_data.to_crs(\"ESRI:54030\")\n", @@ -598,7 +668,7 @@ " ax=ax_bot,\n", " cmap=\"seismic\",\n", " norm=norm_impact,\n", - " edgecolor='black',\n", + " edgecolor=\"black\",\n", " linewidth=0.03,\n", " )\n", " ax_bot.set_axis_off()\n", @@ -610,27 +680,48 @@ " if region_filter:\n", " era5_local_month = era5_local.sel(time=era5_local.time.dt.month == month)\n", "\n", - "\n", - " ref_vals = betas_mmt['mmt'].sel(region=region_filter).sel(age_cohort='age65plus') + 273.15\n", - " betas = betas_mmt['beta_hotonly'] if hotonly else betas_mmt['beta']\n", - " da_temp_bins = betas.sel(region=region_filter).sel(age_cohort='age65plus')\n", - " da_temp_bins['tas_bin'] = da_temp_bins['tas_bin'] + 273.15\n", + " ref_vals = (\n", + " betas_mmt[\"mmt\"].sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + " + 273.15\n", + " )\n", + " betas = betas_mmt[\"beta_hotonly\"] if hotonly else betas_mmt[\"beta\"]\n", + " da_temp_bins = betas.sel(region=region_filter).sel(age_cohort=\"age65plus\")\n", + " da_temp_bins[\"tas_bin\"] = da_temp_bins[\"tas_bin\"] + 273.15\n", "\n", " # Build interpolator for response curve\n", " response_interp = interp1d(\n", - " da_temp_bins.coords['tas_bin'].values,\n", + " da_temp_bins.coords[\"tas_bin\"].values,\n", " da_temp_bins.values,\n", - " bounds_error=False, fill_value=0\n", + " bounds_error=False,\n", + " fill_value=0,\n", " )\n", "\n", - "\n", - "\n", - " ax_hist.hist(bias_adj_iso_all[i].values.flatten(), bins=50, alpha=0.5, label=\"Forecast\", density=density_hist)\n", - " ax_hist.hist(era5_local_month.values.flatten(), bins=50, alpha=0.5, label=\"ERA5 (1995–2025)\", density=density_hist)\n", - " ax_hist.axvline(ref_vals.values, color='black', linewidth=0.5, alpha=0.7, label=\"MMT\")\n", + " ax_hist.hist(\n", + " bias_adj_iso_all[i].values.flatten(),\n", + " bins=50,\n", + " alpha=0.5,\n", + " label=\"Forecast\",\n", + " density=density_hist,\n", + " )\n", + " ax_hist.hist(\n", + " era5_local_month.values.flatten(),\n", + " bins=50,\n", + " alpha=0.5,\n", + " label=\"ERA5 (1995–2025)\",\n", + " density=density_hist,\n", + " )\n", + " ax_hist.axvline(\n", + " ref_vals.values, color=\"black\", linewidth=0.5, alpha=0.7, label=\"MMT\"\n", + " )\n", "\n", " ax_hist2 = ax_hist.twinx()\n", - " ax_hist2.plot(da_temp_bins.coords['tas_bin'], da_temp_bins.values, color='black', linewidth=1.5, label=\"Response curve\")\n", + " ax_hist2.plot(\n", + " da_temp_bins.coords[\"tas_bin\"],\n", + " da_temp_bins.values,\n", + " color=\"black\",\n", + " linewidth=1.5,\n", + " label=\"Response curve\",\n", + " )\n", " ax_hist2.set_ylabel(\"Response\")\n", " ax_hist2.set_ylim(-10, 30)\n", "\n", @@ -644,7 +735,6 @@ " ax3.plot(fc_centers, fc_cdf, label=\"Forecast weighted CDF\")\n", " ax3.plot(era5_centers, era5_cdf, label=\"ERA5 weighted CDF\")\n", "\n", - "\n", " ax_hist.legend(fontsize=7)\n", " ax_hist.set_xlabel(\"Tas (K)\")\n", " ax_hist.set_ylabel(\"Density\")\n", @@ -660,8 +750,8 @@ "fig.canvas.draw()\n", "\n", "pos_top = ax_oct_top.get_position()\n", - "cax_temp = fig.add_axes([pos_top.x1 + 0.01, pos_top.y0+0.02, 0.01, pos_top.height])\n", - "sm_temp = plt.cm.ScalarMappable(cmap='RdBu_r', norm=plt.Normalize(vmin=-5, vmax=5))\n", + "cax_temp = fig.add_axes([pos_top.x1 + 0.01, pos_top.y0 + 0.02, 0.01, pos_top.height])\n", + "sm_temp = plt.cm.ScalarMappable(cmap=\"RdBu_r\", norm=plt.Normalize(vmin=-5, vmax=5))\n", "fig.colorbar(sm_temp, cax=cax_temp, label=\"ΔTas (K)\\nBaseline: 1995-2025\")\n", "\n", "pos_bot = ax_oct_bot.get_position()\n", @@ -677,15 +767,19 @@ "\n", "title_region = f\" {region_filter}\" if region_filter else \"\"\n", "ag_label = \"Age-Weighted\" if age_weight else \"Age Cohort: 65+\"\n", - "label = 'Death Rate' if rate else \"Deaths\"\n", - "label_h = 'Hotonly' if hotonly else 'Net'\n", + "label = \"Death Rate\" if rate else \"Deaths\"\n", + "label_h = \"Hotonly\" if hotonly else \"Net\"\n", "\n", "fig.suptitle(f\"{ag_label} {label_h} {label}:{title_region}\", fontsize=14)\n", "\n", "\n", - "label_f = 'death_rate' if rate else \"deaths\"\n", + "label_f = \"death_rate\" if rate else \"deaths\"\n", "\n", - "fig.savefig(f'{label_h}_{label_f}_agewtd{f'_{region_filter}' if region_filter else ''}.png', dpi=600, bbox_inches=\"tight\")" + "fig.savefig(\n", + " f\"{label_h}_{label_f}_agewtd{f'_{region_filter}' if region_filter else ''}.png\",\n", + " dpi=600,\n", + " bbox_inches=\"tight\",\n", + ")" ] } ],