country scc code

src/dscim/menu/main_recipe.py

@@ -106,6 +106,7 @@
         extrap_formula=None,
         fair_dims=None,
         save_files=None,
+        geography=None,
         **kwargs,
     ):
         if scc_quantiles is None:
@@ -194,6 +195,10 @@
                 "global_consumption_no_pulse",
             ]
 
+        if "country_ISOs" in kwargs:
+            self.countries_mapping = pd.read_csv(kwargs["country_ISOs"])
+            self.countries = self.countries_mapping.MatchedISO.dropna().unique()
+
         super().__init__(
             sector_path=sector_path,
             save_path=save_path,
@@ -203,6 +208,8 @@
             eta=eta,
             subset_dict=subset_dict,
             ce_path=ce_path,
+            geography=geography,
+            kwargs=kwargs,
         )
 
         self.rho = rho
@@ -232,6 +239,7 @@
         self.fair_dims = fair_dims
         self.save_files = save_files
         self.__dict__.update(**kwargs)
+        self.geography = geography
         self.kwargs = kwargs
 
         self.logger = logging.getLogger(__name__)
@@ -477,7 +485,7 @@
 
     @cachedproperty
     @save(name="damage_function_points")
-    def damage_function_points(self) -> pd.DataFrame:
+    def damage_function_points(self, country=None) -> pd.DataFrame:
         """Global damages by RCP/GCM or SLR
 
         Returns
@@ -508,6 +516,37 @@
 
         return df
 
+    def country_damage_function_points(self, country) -> pd.DataFrame:
+        """Global damages by RCP/GCM or SLR
+
+        Returns
+        --------
+            pd.DataFrame
+        """
+        df = self.country_damages_calculation(country)
+
+        if "slr" in df.columns:
+            df = df.merge(self.climate.gmsl, on=["year", "slr"])
+        if "gcm" in df.columns:
+            df = df.merge(self.climate.gmst, on=["year", "gcm", "rcp"])
+
+        # removing illegal combinations from estimation
+        if any([i in df.ssp.unique() for i in ["SSP1", "SSP5"]]):
+            self.logger.info("Dropping illegal model combinations.")
+            for var in [i for i in df.columns if i in ["anomaly", "gmsl"]]:
+                df.loc[
+                    ((df.ssp == "SSP1") & (df.rcp == "rcp85"))
+                    | ((df.ssp == "SSP5") & (df.rcp == "rcp45")),
+                    var,
+                ] = np.nan
+
+        # agriculture lacks ACCESS0-1/rcp85 combo
+        if "agriculture" in self.sector:
+            self.logger.info("Dropping illegal model combinations for agriculture.")
+            df.loc[(df.gcm == "ACCESS1-0") & (df.rcp == "rcp85"), "anomaly"] = np.nan
+
+        return df
+
     def damage_function_calculation(self, damage_function_points, global_consumption):
         """The damage function model fit may be : (1) ssp specific, (2) ssp-model specific, (3) unique across ssp-model.
         This depends on the type of discounting. In each case the input data passed to the fitting functions and the formatting of the returned
@@ -566,47 +605,54 @@
         elif (self.discounting_type == "constant") or (
             "ramsey" in self.discounting_type
         ):
-            for ssp, model in list(
+            loop = list(
                 product(
                     damage_function_points.ssp.unique(),
                     damage_function_points.model.unique(),
                 )
-            ):
-                # Subset dataframe to specific SSP-IAM combination.
-                fit_subset = damage_function_points[
-                    (damage_function_points["ssp"] == ssp)
-                    & (damage_function_points["model"] == model)
-                ]
+            )
+            for country in self.countries if self.geography == "country" else [0]:
+                if self.geography == "country":
+                    damage_function_points = self.country_damage_function_points(
+                        country
+                    )
+                for ssp, model in loop:
+                    # Subset dataframe to specific SSP-IAM combination.
 
-                global_c_subset = global_consumption.sel({"ssp": ssp, "model": model})
+                    selection = (damage_function_points["ssp"] == ssp) & (
+                        damage_function_points["model"] == model
+                    )
 
-                # Fit damage function curves using the data subset
-                damage_function = model_outputs(
-                    damage_function=fit_subset,
-                    formula=self.formula,
-                    type_estimation=self.fit_type,
-                    global_c=global_c_subset,
-                    extrapolation_type=self.ext_method,
-                    quantiles=self.quantreg_quantiles,
-                    year_range=yrs,
-                    year_start_pred=self.ext_subset_end_year + 1,
-                )
+                    fit_subset = damage_function_points[selection]
 
-                # Add variables
-                params = damage_function["parameters"].expand_dims(
-                    dict(
-                        discount_type=[self.discounting_type], ssp=[ssp], model=[model]
-                    )
-                )
+                    subset = dict(ssp=[ssp], model=[model])
 
-                preds = damage_function["preds"].expand_dims(
-                    dict(
-                        discount_type=[self.discounting_type], ssp=[ssp], model=[model]
+                    if self.geography == "country":
+                        subset.update(dict(country=[country]))
+
+                    global_c_subset = global_consumption.sel(subset)
+
+                    # Fit damage function curves using the data subset
+                    damage_function = model_outputs(
+                        damage_function=fit_subset,
+                        formula=self.formula,
+                        type_estimation=self.fit_type,
+                        global_c=global_c_subset,
+                        extrapolation_type=self.ext_method,
+                        quantiles=self.quantreg_quantiles,
+                        year_range=yrs,
+                        year_start_pred=self.ext_subset_end_year + 1,
                     )
-                )
 
-                params_list.append(params)
-                preds_list.append(preds)
+                    subset.update(dict(discount_type=[self.discounting_type]))
+
+                    # Add variables
+                    params = damage_function["parameters"].expand_dims(subset)
+
+                    preds = damage_function["preds"].expand_dims(subset)
+
+                    params_list.append(params)
+                    preds_list.append(preds)
 
         elif "gwr" in self.discounting_type:
             # Fit damage function across all SSP-IAM combinations, as expected
@@ -797,9 +843,15 @@
         """
 
         # Calculate global consumption per capita
-        array_pc = self.global_consumption_calculation(
-            disc_type
-        ) / self.collapsed_pop.sum("region")
+        if self.geography == "country":
+            array_pc = (
+                self.global_consumption_calculation(disc_type)
+                / self.country_econ_vars.pop
+            )
+        else:
+            array_pc = self.global_consumption_calculation(
+                disc_type
+            ) / self.collapsed_pop.sum("region")
 
         if self.NAME == "equity":
             # equity recipe's growth is capped to
@@ -840,7 +892,10 @@
 
             # holding population constant
             # from 2100 to 2300 with 2099 values
-            pop = self.collapsed_pop.sum("region")
+            if self.geography == "country":
+                pop = self.country_econ_vars.pop
+            else:
+                pop = self.collapsed_pop.sum("region")
             pop = pop.reindex(
                 year=range(pop.year.min().values, self.ext_end_year + 1),
                 method="ffill",

src/dscim/menu/risk_aversion.py

@@ -1,6 +1,7 @@
 import pandas as pd
 import xarray as xr
 from dscim.menu.main_recipe import MainRecipe
+from dscim.descriptors import cachedproperty
 
 
 class RiskAversionRecipe(MainRecipe):
@@ -38,7 +39,7 @@
 
         return ce_array
 
-    @property
+    @cachedproperty
     def calculated_damages(self) -> xr.DataArray:
         """Calculate damages (difference between CEs)"""
         return self.ce_no_cc - self.ce_cc
@@ -62,6 +63,43 @@
 
         return df
 
+    def country_damages_calculation(self, country) -> pd.DataFrame:
+        """Aggregate damages to country level
+
+        Returns
+        --------
+            pd.DataFrame
+        """
+
+        self.logger.info(f"Calculating damages for {country}")
+        dams_collapse = self.calculated_damages * self.collapsed_pop
+
+        territories = []
+        for ii, row in self.countries_mapping.iterrows():
+            if row["MatchedISO"] == country:
+                territories.append(row["ISO"])
+        self.logger.debug(f"Including territories {territories}")
+
+        dams_collapse = dams_collapse.sel(
+            region=[
+                i
+                for i in dams_collapse.region.values
+                if any(j in i for j in territories)
+            ]
+        )
+
+        dams_collapse = (
+            dams_collapse.sum(dim="region").to_dataframe("damages").reset_index()
+        )
+
+        if "gwr" in self.discounting_type:
+            dams_collapse = dams_collapse.assign(
+                ssp=str(list(self.gdp.ssp.values)),
+                model=str(list(self.gdp.model.values)),
+            )
+
+        return dams_collapse
+
     def global_consumption_calculation(self, disc_type):
         """Calculate global consumption
 
@@ -71,7 +109,10 @@
         """
 
         if (disc_type == "constant") or ("ramsey" in disc_type):
-            global_cons_no_cc = self.gdp.sum(dim=["region"])
+            if self.geography == "country":
+                global_cons_no_cc = self.country_econ_vars.gdp
+            else:
+                global_cons_no_cc = self.gdp.sum(dim=["region"])
 
         elif disc_type == "constant_model_collapsed":
             global_cons_no_cc = self.gdp.sum(dim=["region"]).mean(dim=["model"])

src/dscim/menu/simple_storage.py

@@ -306,7 +306,12 @@
         histclim=None,
         ce_path=None,
         subset_dict=None,
+        geography=None,
+        **kwargs,
     ):
+        if geography is None:
+            geography = "global"
+
         self.sector_path = sector_path
         self.save_path = save_path
         self.gdppc_bottom_code = gdppc_bottom_code
@@ -317,6 +322,13 @@
         self.histclim = histclim
         self.ce_path = ce_path
         self.eta = eta
+        self.__dict__.update(**kwargs)
+        self.geography = geography
+        self.kwargs = kwargs
+
+        if "country_ISOs" in kwargs:
+            self.countries_mapping = pd.read_csv(kwargs["country_ISOs"])
+            self.countries = self.countries_mapping.MatchedISO.dropna().unique()
 
         self.logger = logging.getLogger(__name__)
 
@@ -368,6 +380,32 @@
             )
         return raw
 
+    @property
+    def country_econ_vars(self):
+        mapping_dict = {}
+        for ii, row in self.countries_mapping.iterrows():
+            mapping_dict[row["ISO"]] = row["MatchedISO"]
+            if row["MatchedISO"] == "nan":
+                mapping_dict[row["ISO"]] = "nopop"
+
+        socioec = self.econ_vars.econ_vars
+        socioec = socioec.assign_coords(
+            {"region": [region[:3] for region in socioec.region.values]}
+        )
+        new_ISOs = []
+        for ISO in socioec.region.values:
+            new_ISOs.append(mapping_dict[ISO])
+        raw = (
+            socioec.assign_coords({"region": new_ISOs})
+            .groupby("region")
+            .sum()
+            .rename({"region": "country"})
+            .drop_sel(country="nan")
+        )
+        raw = self.cut(raw, end_year=2099)
+
+        return raw
+
     @cachedproperty
     def gdp(self):
         return self.cut_econ_vars.gdp

src/dscim/utils/rff.py

@@ -364,15 +364,19 @@
     weights_path,
     pulse_year,
     mask,
+    geography=None,
 ):
     """Wrapper function for `weight_df()`."""
 
     # ssp GDP for fractionalizing damage functions
-    ssp_gdp = xr.open_zarr(ssp_gdp, consolidated=True).sum("region").gdp
-
-    # get RFF data
-    region = "USA" if USA else "world"
-    rff_gdp = xr.open_dataset(rff_gdp).sel(region=region, drop=True).gdp
+    if geography == "country":
+        ssp_gdp = xr.open_zarr(ssp_gdp, consolidated=True).gdp
+        rff_gdp = xr.open_dataset(rff_gdp).gdp
+    else:
+        ssp_gdp = xr.open_zarr(ssp_gdp, consolidated=True).sum("region").gdp
+        # get RFF data
+        region = "USA" if USA else "world"
+        rff_gdp = xr.open_dataset(rff_gdp).sel(region=region, drop=True).gdp
 
     # get global consumption factors to extrapolate damage function
     factors = rff_gdp.sel(year=slice(2100, 2300)) / rff_gdp.sel(year=2099)