diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e89c900..6f3b21c9 100644 --- a/CHANGELOG.md +++ b/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)) diff --git a/src/dscim/preprocessing/input_damages.py b/src/dscim/preprocessing/input_damages.py index dc36b7db..33af9232 100644 --- a/src/dscim/preprocessing/input_damages.py +++ b/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, diff --git a/src/dscim/preprocessing/mortality_inputs_archived.py b/src/dscim/preprocessing/mortality_inputs_archived.py deleted file mode 100644 index e3dde78f..00000000 --- a/src/dscim/preprocessing/mortality_inputs_archived.py +++ /dev/null @@ -1,66 +0,0 @@ -import xarray as xr -from dscim.utils.functions import gcms -from pathlib import Path -from dscim.preprocessing.input_damages import prep_mortality_damages - -######################### -# MORTALITY V4 -######################### - -prep_mortality_damages( - gcms=gcms(), - paths=dict( - delta_deaths=[ - Path( - f"/shares/gcp/outputs/mortality/impacts-darwin-montecarlo-damages/vsl_popavg/mortality_damages_IR_batch{i}.nc4" - ) - for i in range(15) - ], - delta_costs=[ - Path( - f"/shares/gcp/outputs/mortality/impacts-darwin-montecarlo-damages/mortality_damages_IR_batch{i}.nc4" - ) - for i in range(15) - ], - histclim_deaths=[ - Path( - f"/shares/gcp/outputs/mortality/impacts-darwin-montecarlo-damages/vsl_popavg/mortality_damages_IR_batch{i}_histclim.nc4" - ) - for i in range(15) - ], - ), - vars=dict( - delta_deaths="monetized_deaths_vsl_epa_popavg", - delta_costs="monetized_costs_vsl_epa_scaled", - histclim_deaths="monetized_deaths_vsl_epa_popavg", - ), - outpath="/shares/gcp/integration/float32/sectoral_ir_damages/mortality_data/impacts-darwin-montecarlo-damages-v4.zarr", -) - -######################### -# MORTALITY V5 -######################### - -# prep_mortality_damages( -# gcms=gcms(), -# paths=dict( -# delta_deaths = [ -# Path(f"/shares/gcp/estimation/mortality/release_2020/data/3_valuation/impact_region/row/mortality_damages_IR_batch{i}.nc4") -# for i in range(15) -# ], -# delta_costs = [ -# Path(f"/shares/gcp/outputs/mortality/impacts-darwin-montecarlo-damages/mortality_damages_IR_batch{i}.nc4") -# for i in range(15) -# ], -# histclim_deaths = [ -# Path(f"/shares/gcp/estimation/mortality/release_2020/data/3_valuation/impact_region/row/mortality_damages_IR_batch{i}.nc4") -# for i in range(15) -# ], -# ), -# vars=dict( -# delta_deaths = 'monetized_deaths_vsl_epa_scaled', -# delta_costs = 'monetized_costs_vsl_epa_scaled', -# histclim_deaths = 'monetized_histclim_deaths_vsl_epa_scaled', -# ), -# outpath="/shares/gcp/integration/float32/sectoral_ir_damages/mortality_data/impacts-darwin-montecarlo-damages-v5.zarr", -# )