Merge pull request #1 from ClimateImpactLab/functions

PR for functions branch

Author: kitschwarz

dscim/diagnostics/damage_inc_share_boxplots.py

@@ -0,0 +1,118 @@
+import xarray as xr
+import pandas as pd
+import os, sys, yaml
+
+USER = os.getenv("USER")
+from dscim.utils.functions import get_model_weights, US_territories
+from impactlab_tools.utils.weighting import weighted_quantile_xr
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def income_boxplot(
+    damages,
+    delta_var,
+    socioec,
+    quantile,
+    sector,
+    rcp,
+    ssp,
+    model,
+    year,
+    USA,
+    output,
+):
+
+    if quantile == "quintile":
+        bin_breaks = np.linspace(0, 1, 6)
+    elif quantile == "decile":
+        bin_breaks = np.linspace(0, 1, 11)
+
+    # open files and select relevant data
+    socioec = xr.open_zarr(socioec).sel(ssp=ssp, model=model)
+    mean_damages = xr.open_zarr(damages).sel(ssp=ssp, model=model, year=year, rcp=rcp)[
+        delta_var
+    ]
+
+    if USA == True:
+        # subset US data
+        US_IRs = [
+            i for i in socioec.region.values if any(j in i for j in US_territories())
+        ]
+        socioec = socioec.sel(region=US_IRs)
+        mean_damages = mean_damages.sel(region=US_IRs)
+
+    if "CAMEL" not in sector:
+        # reduce damages for non-CAMEL sectors (CAMEL is generated differently and is thus pre-collapsed)
+        weights = get_model_weights(rcp)
+        mean_damages = mean_damages.mean("batch").weighted(weights).mean("gcm")
+
+    # get damages as a share of income
+    damages_inc_share = (
+        mean_damages
+        / socioec.gdppc.sel(year=2020).rename("damages_inc_share").to_dataset()
+    )
+    damages_inc_share["present_gdppc"] = socioec.gdppc.sel(year=2020)
+    damages_inc_share["present_pop"] = socioec["pop"].sel(year=2020)
+    damages_inc_share = damages_inc_share.dropna("region")
+
+    # pop-weighted gdppc quantiles (determine the bin cutoffs)
+    bins = weighted_quantile_xr(
+        damages_inc_share.present_gdppc,
+        bin_breaks,
+        damages_inc_share.present_pop,
+        "region",
+    ).values
+
+    groupby = damages_inc_share.groupby_bins("present_gdppc", bins)
+
+    plot_bins = []
+    quantiles = [0.05, 0.25, 0.5, 0.75, 0.95]
+
+    # get quantiles of damages_inc_share for each gdppc bin
+    for i, j in groupby:
+
+        bin_no = [round(i, 1) for i in bins].index(round(i.left, 1))
+
+        values = list(
+            weighted_quantile_xr(
+                j.damages_inc_share, quantiles, j.present_pop, "region"
+            ).values
+        )
+
+        values.append(j.damages_inc_share.weighted(j.present_pop).mean())
+
+        plot_bins.append((bin_no + 1, values))
+
+    # plot them
+    fig, ax = plt.subplots(figsize=(15, 10), sharey=True)
+
+    stats = []
+
+    for bin, s in sorted(plot_bins):
+
+        stats.append(
+            {
+                "label": bin,
+                "mean": s[5],
+                "med": s[2],
+                "q1": s[1],
+                "q3": s[3],
+                "whislo": s[0],
+                "whishi": s[4],
+                "fliers": [],
+            }
+        )
+
+    bxp = ax.bxp(stats, showmeans=True, meanline=False)
+    ax.set_ylabel(f"{year} damages as a share of present-day GDP")
+    ax.set_xlabel("Decile of present-day GDP per capita")
+    ax.spines["top"].set_visible(False)
+    ax.spines["right"].set_visible(False)
+
+    os.makedirs(output, exist_ok=True)
+    plt.savefig(
+        f"{output}/boxplot_{sector}_domestic-{USA}_{rcp}_{ssp}_{model}_{quantile}_{year}.png",
+        bbox_inches="tight",
+        dpi=300,
+    )

dscim/diagnostics/var_timeseries.py

@@ -24,13 +24,13 @@ def get_rff_id(
 ):
 
     if mask in ["epa_mask", "unmasked"]:
-        results = f"{results_root}/{sector}/{pulse_year}/unmasked_None/"
+        pass
     elif mask == "gdppc_mask":
-        results = f"{results_root}/{sector}/{pulse_year}/gdppc_emissions_q0.01_q0.99"
+        print("gdppc_mask deprecated")
 
     sccs = (
         xr.open_dataset(
-            f"{results}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_sccs.nc4"
+            f"{results_root}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_sccs.nc4"
         )
         .sel(
             weitzman_parameter="0.5",
@@ -143,11 +143,8 @@ def rff_timeseries(
     runid_root,
     discrate=0.02,
     pulse_year=2020,
-    results_mask="unmasked_None",
 ):
 
-    results = f"{results_root}/{sector}/{pulse_year}/{results_mask}"
-
     # index
     rff_ids = (
         pd.read_csv(
@@ -163,7 +160,7 @@ def rff_timeseries(
     # sccs
     sccs = (
         xr.open_dataset(
-            f"{results}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_sccs.nc4"
+            f"{results_root}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_sccs.nc4"
         )
         .sel(
             weitzman_parameter="0.5", gas=gas, fair_aggregation="uncollapsed", drop=True
@@ -174,7 +171,7 @@ def rff_timeseries(
     # marginal damages
     damages = (
         xr.open_zarr(
-            f"{results}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_marginal_damages.zarr"
+            f"{results_root}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_marginal_damages.zarr"
         )
         .sel(weitzman_parameter="0.5", gas=gas, drop=True)
         .marginal_damages
@@ -183,7 +180,7 @@ def rff_timeseries(
     # discount factors
     df = (
         xr.open_zarr(
-            f"{results}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_discount_factors.zarr"
+            f"{results_root}/{recipe}_{disc}_eta{eta}_rho{rho}_uncollapsed_discount_factors.zarr"
         )
         .sel(weitzman_parameter="0.5", gas=gas, drop=True)
         .discount_factor.rename("discount_factors")

dscim/menu/equity.py

@@ -17,7 +17,7 @@ def ce_cc_calculation(self) -> xr.DataArray:
         else:
             dims = ["region"]
 
-        ce_cc = self.risk_aversion_damages("ce_cc").ce_cc
+        ce_cc = self.risk_aversion_damages("cc").cc
 
         # reindex to make sure all regions are being calculated
         if len(ce_cc.region.values) < len(self.gdppc.region.values):
@@ -32,15 +32,15 @@ def ce_cc_calculation(self) -> xr.DataArray:
             eta=self.eta,
         )
 
-        return ce_array.rename("ce_cc")
+        return ce_array.rename("cc")
 
     def ce_no_cc_calculation(self) -> xr.DataArray:
         if "gwr" in self.discounting_type:
             dims = ["ssp", "model", "region"]
         else:
             dims = ["region"]
 
-        ce_no_cc = self.risk_aversion_damages("ce_no_cc").ce_no_cc
+        ce_no_cc = self.risk_aversion_damages("no_cc").no_cc
         # reindex to make sure all regions are being calculated
         if len(ce_no_cc.region.values) < len(self.gdppc.region.values):
             ce_no_cc = ce_no_cc.reindex({"region": self.gdppc.region.values})
@@ -54,7 +54,7 @@ def ce_no_cc_calculation(self) -> xr.DataArray:
             eta=self.eta,
         )
 
-        return ce_no_cc_array.rename("ce_no_cc")
+        return ce_no_cc_array.rename("no_cc")
 
     @property
     def calculated_damages(self) -> xr.DataArray:

dscim/menu/main_recipe.py

@@ -71,12 +71,12 @@ class MainRecipe(StackedDamages, ABC):
 
     def __init__(
         self,
-        sector_path,
-        save_path,
         econ_vars,
         climate_vars,
         sector,
         formula,
+        sector_path=None,
+        save_path=None,
         rho=0.00461878399,
         eta=1.421158116,
         fit_type="ols",
@@ -92,7 +92,7 @@ def __init__(
         gdppc_bottom_code=39.39265060424805,
         scc_quantiles=[0.05, 0.17, 0.25, 0.5, 0.75, 0.83, 0.95],
         scenario_dimensions=None,
-        weitzman_parameter=[0.1],
+        weitzman_parameter=[0.1, 0.5],
         fair_aggregation=["ce", "mean", "gwr_mean", "median", "median_params"],
         filename_suffix="",
         discrete_discounting=False,

dscim/menu/risk_aversion.py

@@ -17,7 +17,7 @@ def ce_cc_calculation(self) -> xr.DataArray:
         -------
              xr.DataArray
         """
-        ce_array = self.risk_aversion_damages("ce_cc").ce_cc
+        ce_array = self.risk_aversion_damages("cc").cc
 
         # for GWR options, take the CE over growth models
         if "gwr" in self.discounting_type:
@@ -32,7 +32,7 @@ def ce_no_cc_calculation(self) -> xr.DataArray:
         -------
             xr.DataArray
         """
-        ce_array = self.risk_aversion_damages("ce_no_cc").ce_no_cc
+        ce_array = self.risk_aversion_damages("no_cc").no_cc
 
         if "gwr" in self.discounting_type:
             ce_array = self.ce(ce_array, dims=["ssp", "model"])

dscim/menu/simple_storage.py

@@ -46,10 +46,10 @@ def __init__(
         gmsl_path,
         gmst_fair_path,
         damages_pulse_conversion_path,
+        pulse_year,
         gmsl_fair_path=None,
         ecs_mask_path=None,
         ecs_mask_name=None,
-        pulse_year=2020,
         base_period=[2001, 2010],
         emission_scenarios=["ssp119", "ssp126", "ssp245", "ssp460", "ssp370", "ssp585"],
         gases=["CO2_Fossil", "CH4", "N2O"],
@@ -388,16 +388,15 @@ def gdppc(self):
     def adding_up_damages(self):
         """This property calls pre-calculated adding-up IR-level 'mean' over batches."""
 
-        mean_cc = f"{self.ce_path}/adding_up_ce_cc_eta{self.eta}.zarr"
-        mean_no_cc = f"{self.ce_path}/adding_up_ce_no_cc_eta{self.eta}.zarr"
+        mean_cc = f"{self.ce_path}/adding_up_cc.zarr"
+        mean_no_cc = f"{self.ce_path}/adding_up_no_cc.zarr"
 
         if os.path.exists(mean_cc) and os.path.exists(mean_no_cc):
             self.logger.info(
                 f"Adding up aggregated damages found at {mean_cc}, {mean_no_cc}. These are being loaded..."
             )
             damages = (
-                (xr.open_zarr(mean_no_cc).ce_no_cc - xr.open_zarr(mean_cc).ce_cc)
-                * self.pop
+                (xr.open_zarr(mean_no_cc).no_cc - xr.open_zarr(mean_cc).cc) * self.pop
             ).sum("region")
         else:
             raise NotImplementedError(
@@ -410,7 +409,7 @@ def risk_aversion_damages(self, ce_type):
 
         Parameters
         ----------
-        ce_type : either `ce_no_cc` or `ce_cc`
+        ce_type : either `no_cc` or `cc`
 
         Returns
         -------