Merge pull request #74 from ClimateImpactLab/new_mortality_format

Change mortality input file to be compatible with new format

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - Add functions used for producing RFF weights. ([PR #66](https://github.com/ClimateImpactLab/dscim/pull/66), [@davidrzhdu](https://github.com/davidrzhdu))
 ### Changed
+- Changed `prep_mortality_damages` function to work with new format mortality outputs. ([PR #74](https://github.com/ClimateImpactLab/dscim/pull/65), [@JMGilbert](https://github.com/JMGilbert))
 - Included US territories in damages and economic variable subsetting. ([PR #78](https://github.com/ClimateImpactLab/dscim/pull/78), [@JMGilbert](https://github.com/JMGilbert))
 - Changed format of `eta_rhos` to allow for multiple values of `rho` for the same `eta`. ([PR #65](https://github.com/ClimateImpactLab/dscim/pull/65), [@JMGilbert](https://github.com/JMGilbert))
 - Remove diagnostics module. ([PR #60](https://github.com/ClimateImpactLab/dscim/pull/60), [@JMGilbert](https://github.com/JMGilbert))

src/dscim/preprocessing/input_damages.py

@@ -656,36 +656,75 @@ def prep_mortality_damages(
     paths,
     vars,
     outpath,
+    mortality_version,
     path_econ="/shares/gcp/integration/float32/dscim_input_data/econvars/zarrs/integration-econ-bc39.zarr",
 ):
-    print(
-        "This function only works on mortality_v4 and mortality_v5 damages from the mortality repo's valuation. Earlier versions of mortality contain different variable definitions (per capita, not per capita, with or without histclim subtracted off."
-    )
-
     ec = EconVars(path_econ=path_econ)
 
     # longest-string gcm has to be processed first so the coordinate is the right str length
     gcms = sorted(gcms, key=len, reverse=True)
 
-    for i, gcm in enumerate(gcms):
-        print(gcm, i, "/", len(gcms))
+    if mortality_version == 0:
+        scaling_deaths = "epa_scaled"
+        scaling_costs = "epa_scaled"
+        valuation = "vly"
+    elif mortality_version == 1:
+        scaling_deaths = "epa_iso_scaled"
+        scaling_costs = "epa_scaled"
+        valuation = "vsl"
+    elif mortality_version == 4:
+        scaling_deaths = "epa_popavg"
+        scaling_costs = "epa_scaled"
+        valuation = "vsl"
+    elif mortality_version == 5:
+        scaling_deaths = "epa_row"
+        scaling_costs = "epa_scaled"
+        valuation = "vsl"
+    else:
+        raise ValueError("Mortality version not valid: ", str(mortality_version))
+
+    # We set 0 population to infinity so that per capita damages are 0 in locations with 0 population
+    pop = ec.econ_vars.pop.load()
+    pop = xr.where(pop == 0, np.Inf, pop)
 
-        def prep(ds, gcm=gcm):
-            return ds.sel(gcm=gcm).drop("gcm")
+    for i, gcm in enumerate(gcms):
+        print(gcm, i + 1, "/", len(gcms))
 
         data = {}
         for var, name in vars.items():
-            data[var] = xr.open_mfdataset(paths[var], preprocess=prep, parallel=True)[
-                name
-            ]
+
+            def prep(
+                ds,
+                gcm=gcm,
+                scaling_deaths=scaling_deaths,
+                scaling_costs=scaling_costs,
+                valuation=valuation,
+                name=name,
+            ):
+                return ds.sel(
+                    gcm=gcm,
+                    scaling=[scaling_deaths, scaling_costs],
+                    valuation=valuation,
+                ).drop(["gcm", "valuation"])
+
+            data = xr.open_mfdataset(
+                paths, preprocess=prep, parallel=True, engine="zarr"
+            )
 
         damages = xr.Dataset(
             {
                 "delta": (
-                    data["delta_deaths"] - data["histclim_deaths"] + data["delta_costs"]
+                    data[vars["delta_deaths"]].sel(scaling=scaling_deaths, drop=True)
+                    - data[vars["histclim_deaths"]].sel(
+                        scaling=scaling_deaths, drop=True
+                    )
+                    + data[vars["delta_costs"]].sel(scaling=scaling_costs, drop=True)
+                )
+                / pop,
+                "histclim": data[vars["histclim_deaths"]].sel(
+                    scaling=scaling_deaths, drop=True
                 )
-                / ec.econ_vars.pop.load(),
-                "histclim": data["histclim_deaths"] / ec.econ_vars.pop.load(),
+                / pop,
             }
         ).expand_dims({"gcm": [gcm]})
 
@@ -697,6 +736,14 @@ def prep(ds, gcm=gcm):
         # convert to EPA VSL
         damages = damages * 0.90681089
 
+        for v in list(damages.coords.keys()):
+            if damages.coords[v].dtype == object:
+                damages.coords[v] = damages.coords[v].astype("unicode")
+
+        for v in list(damages.variables.keys()):
+            if damages[v].dtype == object:
+                damages[v] = damages[v].astype("unicode")
+
         if i == 0:
             damages.to_zarr(
                 outpath,