From dc0389b183291d59ed937fe229407cf4f7c19f05 Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Sun, 5 Jan 2025 07:22:39 -0800 Subject: [PATCH 01/10] include floating point adders to lslr dim in surge lookup --- pyCIAM/surge/lookup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyCIAM/surge/lookup.py b/pyCIAM/surge/lookup.py index 43bb7bb..935e626 100644 --- a/pyCIAM/surge/lookup.py +++ b/pyCIAM/surge/lookup.py @@ -130,8 +130,8 @@ def _get_lslr_rhdiff_range( mins.append(min_lslr) maxs.append(max_lslr) - min_lslr = xr.concat(mins, dim="tmp").min("tmp") - max_lslr = xr.concat(maxs, dim="tmp").max("tmp") + min_lslr = xr.concat(mins, dim="tmp").min("tmp") - 2 * np.finfo("float32").eps + max_lslr = xr.concat(maxs, dim="tmp").max("tmp") + 2 * np.finfo("float32").eps at = _get_planning_period_map(lslr.year, at_start) From c2d31e7cffd7262e90c9e4b0e15939fb76898985 Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Sun, 5 Jan 2025 07:23:22 -0800 Subject: [PATCH 02/10] do better job checking that surge lookup interp only results in nans for expected low lslr --- pyCIAM/run.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index 01aee3a..7df2c60 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -341,21 +341,27 @@ def calc_costs( ) if this_surge_lookup.sum() == 0: continue - surge_noadapt.append( + this_surge_noadapt = ( this_surge_lookup.sel(adapttype="retreat", drop=True) .interp( lslr=lslr.sel(seg=seg), rh_diff=rh_diff_noadapt.sel(seg=seg), assume_sorted=True, - kwargs={"fill_value": 0}, ) .reset_coords(drop=True) .expand_dims(seg=[seg]) ) + # ensure nans are only at the beginning + assert this_surge_noadapt.isnull().sum( + "lslr" + ) == this_surge_noadapt.notnull().argmax("lslr") + + surge_noadapt.append(this_surge_noadapt.fillna(0)) + surge_adapt = [] for adapttype in this_surge_lookup.adapttype.values: - surge_adapt.append( + this_surge_adapt = ( this_surge_lookup.sel(adapttype=adapttype) .interp( lslr=lslr.sel(seg=seg), @@ -363,10 +369,16 @@ def calc_costs( adapttype=adapttype, seg=seg, drop=True ), assume_sorted=True, - kwargs={"fill_value": 0}, ) .reset_coords(drop=True) ) + + # ensure nans are only at the beginning + assert this_surge_adapt.isnull().sum( + "lslr" + ) == this_surge_adapt.notnull().argmax("lslr") + + surge_adapt.append(this_surge_adapt.fillna(0)) surge.append( xr.concat(surge_adapt, dim=this_surge_lookup.adapttype).expand_dims( seg=[seg] From 54433e06bd10fa0f958ca7401ffe26342882d82e Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Mon, 6 Jan 2025 04:33:34 +0000 Subject: [PATCH 03/10] fix check for correctly interpolated surge impacts --- pyCIAM/run.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index 7df2c60..b3c5d69 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -339,6 +339,7 @@ def calc_costs( .reset_coords(drop=True) .frac_losses.rename(lslr_by_seg="lslr", rh_diff_by_seg="rh_diff") ) + lslr_too_low = lslr.sel(seg=seg) < this_surge_lookup.lslr.min() if this_surge_lookup.sum() == 0: continue this_surge_noadapt = ( @@ -353,9 +354,7 @@ def calc_costs( ) # ensure nans are only at the beginning - assert this_surge_noadapt.isnull().sum( - "lslr" - ) == this_surge_noadapt.notnull().argmax("lslr") + assert (this_surge_noadapt.notnull() | lslr_too_low).all() surge_noadapt.append(this_surge_noadapt.fillna(0)) @@ -374,9 +373,7 @@ def calc_costs( ) # ensure nans are only at the beginning - assert this_surge_adapt.isnull().sum( - "lslr" - ) == this_surge_adapt.notnull().argmax("lslr") + assert (this_surge_adapt.notnull() | lslr_too_low).all() surge_adapt.append(this_surge_adapt.fillna(0)) surge.append( @@ -1078,7 +1075,6 @@ def execute_pyciam( # determine whether to check for finished jobs if output_path is None: check = False - tmp_output_path = None else: check = True From 7241da86c435827dc847a95c1c2804ae720bcb46 Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Sun, 5 Jan 2025 20:34:32 -0800 Subject: [PATCH 04/10] improve comments in _load_lslr_for_ciam --- pyCIAM/io.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyCIAM/io.py b/pyCIAM/io.py index f535864..7a0d196 100644 --- a/pyCIAM/io.py +++ b/pyCIAM/io.py @@ -310,14 +310,16 @@ def _load_lslr_for_ciam( .set_index(scen_mc=ix_names) ) - # interpolate to yearly + # add on base year where slr is 0 slr_out = slr_out.reindex( year=np.concatenate(([slr_0_year], slr.year.values)), fill_value=0, ) + # interpolate to desired years if interp_years is not None: slr_out = slr_out.interp(year=interp_years) + return slr_out From d2f07d42726393e3ee41669b832ee26d79f2cb9b Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Tue, 7 Jan 2025 01:25:00 +0000 Subject: [PATCH 05/10] more informative surge interpolation errors, and catch failed interpolation --- pyCIAM/run.py | 45 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 36 insertions(+), 9 deletions(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index b3c5d69..a4c37ec 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -332,6 +332,24 @@ def calc_costs( # interpolation with 0's surge_noadapt = [] surge = [] + + def _check_vals(rh_diff_arr, lslr_arr, seg): + error_str = ( + "{0} value is {1} than {2} in storm damage lookup " + f"table for segment {seg}. Please investigate why this is. You may " + "need to re-generate the surge lookup table (or perhaps the `refA` " + "(initial adaptation) table if it was generated for a different " + "set of SLR scenarios or years." + ) + if rh_diff_arr.max() > this_surge_lookup.rh_diff.max(): + raise ValueError(error_str.format("rh_diff", "higher", "maximum")) + if rh_diff_arr.min() < this_surge_lookup.rh_diff.min(): + raise ValueError(error_str.format("rh_diff", "lower", "minimum")) + if lslr_arr.max() > this_surge_lookup.lslr.max(): + raise ValueError(error_str.format("lslr", "higher", "maximum")) + if lslr_arr.min() < this_surge_lookup.lslr.min(): + raise ValueError(error_str.format("lslr", "lower", "minimum")) + for seg in inputs.seg.values: this_surge_lookup = ( surge_lookup.sel(seg=seg) @@ -342,11 +360,16 @@ def calc_costs( lslr_too_low = lslr.sel(seg=seg) < this_surge_lookup.lslr.min() if this_surge_lookup.sum() == 0: continue + + this_rh_diff_noadapt = rh_diff_noadapt.sel(seg=seg, drop=True) + this_lslr = lslr.sel(seg=seg, drop=True) + _check_vals(this_rh_diff_noadapt, this_lslr, seg) + this_surge_noadapt = ( this_surge_lookup.sel(adapttype="retreat", drop=True) .interp( - lslr=lslr.sel(seg=seg), - rh_diff=rh_diff_noadapt.sel(seg=seg), + lslr=this_lslr, + rh_diff=this_rh_diff_noadapt, assume_sorted=True, ) .reset_coords(drop=True) @@ -354,18 +377,22 @@ def calc_costs( ) # ensure nans are only at the beginning - assert (this_surge_noadapt.notnull() | lslr_too_low).all() + assert (this_surge_noadapt.notnull() | lslr_too_low).all(), seg surge_noadapt.append(this_surge_noadapt.fillna(0)) surge_adapt = [] + + this_rh_diff_adapt = rh_diff.sel(seg=seg, drop=True) + _check_vals(this_rh_diff_noadapt, this_lslr, seg) + for adapttype in this_surge_lookup.adapttype.values: this_surge_adapt = ( this_surge_lookup.sel(adapttype=adapttype) .interp( - lslr=lslr.sel(seg=seg), - rh_diff=rh_diff.sel( - adapttype=adapttype, seg=seg, drop=True + lslr=this_lslr, + rh_diff=this_rh_diff_adapt.sel( + adapttype=adapttype, drop=True ), assume_sorted=True, ) @@ -373,7 +400,7 @@ def calc_costs( ) # ensure nans are only at the beginning - assert (this_surge_adapt.notnull() | lslr_too_low).all() + assert (this_surge_adapt.notnull() | lslr_too_low).all(), seg surge_adapt.append(this_surge_adapt.fillna(0)) surge.append( @@ -1157,7 +1184,7 @@ def execute_pyciam( storage_options=storage_options, ) # block on this calculation - wait(surge_futs) + client.gather(surge_futs) ############################### # define temporary output store @@ -1380,7 +1407,7 @@ def execute_pyciam( ############################### # Rechunk and save final ############################### - wait(ciam_futs_2.tolist()) + client.gather(ciam_futs_2.tolist()) assert [f.status == "finished" for f in ciam_futs_2.tolist()] client.cancel(ciam_futs_2) del ciam_futs_2 From 5f0c2c4c6f05cdf3e624fb631976ed800c006a79 Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Tue, 7 Jan 2025 01:39:38 +0000 Subject: [PATCH 06/10] fix error catching in surge lookup --- pyCIAM/run.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index a4c37ec..e4e395c 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -341,14 +341,14 @@ def _check_vals(rh_diff_arr, lslr_arr, seg): "(initial adaptation) table if it was generated for a different " "set of SLR scenarios or years." ) + # lslr is allowed to be below min surge lookup value b/c we created + # lookup such that < min value will be 0 impact if rh_diff_arr.max() > this_surge_lookup.rh_diff.max(): raise ValueError(error_str.format("rh_diff", "higher", "maximum")) if rh_diff_arr.min() < this_surge_lookup.rh_diff.min(): raise ValueError(error_str.format("rh_diff", "lower", "minimum")) if lslr_arr.max() > this_surge_lookup.lslr.max(): raise ValueError(error_str.format("lslr", "higher", "maximum")) - if lslr_arr.min() < this_surge_lookup.lslr.min(): - raise ValueError(error_str.format("lslr", "lower", "minimum")) for seg in inputs.seg.values: this_surge_lookup = ( From fdffaf098b993014a39969c78b081719c618dfcd Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Tue, 7 Jan 2025 16:09:56 +0000 Subject: [PATCH 07/10] allow high rh_diff for surge_noadapt lookup --- pyCIAM/run.py | 37 ++++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index e4e395c..ef5018c 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -333,7 +333,13 @@ def calc_costs( surge_noadapt = [] surge = [] - def _check_vals(rh_diff_arr, lslr_arr, seg): + def _check_vals(rh_diff_arr, lslr_arr, seg, check_rh_diff_max=True): + # for surge_noadapt, rh_diff can be greater than max. This happens when + # a segment chooses retreat10000 in their refA (initial adaptation) and + # also has a declining no-climate-change SLR scenario (e.g. for + # no-climate-change scenarios). However, we automatically assign 0 + # expected storm impacts for retreat10000, so this is acceptable and + # does not need to exist in the surge lookup table. error_str = ( "{0} value is {1} than {2} in storm damage lookup " f"table for segment {seg}. Please investigate why this is. You may " @@ -343,7 +349,10 @@ def _check_vals(rh_diff_arr, lslr_arr, seg): ) # lslr is allowed to be below min surge lookup value b/c we created # lookup such that < min value will be 0 impact - if rh_diff_arr.max() > this_surge_lookup.rh_diff.max(): + if ( + check_rh_diff_max + and rh_diff_arr.max() > this_surge_lookup.rh_diff.max() + ): raise ValueError(error_str.format("rh_diff", "higher", "maximum")) if rh_diff_arr.min() < this_surge_lookup.rh_diff.min(): raise ValueError(error_str.format("rh_diff", "lower", "minimum")) @@ -357,13 +366,16 @@ def _check_vals(rh_diff_arr, lslr_arr, seg): .reset_coords(drop=True) .frac_losses.rename(lslr_by_seg="lslr", rh_diff_by_seg="rh_diff") ) - lslr_too_low = lslr.sel(seg=seg) < this_surge_lookup.lslr.min() if this_surge_lookup.sum() == 0: continue this_rh_diff_noadapt = rh_diff_noadapt.sel(seg=seg, drop=True) this_lslr = lslr.sel(seg=seg, drop=True) - _check_vals(this_rh_diff_noadapt, this_lslr, seg) + _check_vals( + this_rh_diff_noadapt, this_lslr, seg, check_rh_diff_max=False + ) + + lslr_too_low = this_lslr < this_surge_lookup.lslr.min() this_surge_noadapt = ( this_surge_lookup.sel(adapttype="retreat", drop=True) @@ -376,8 +388,13 @@ def _check_vals(rh_diff_arr, lslr_arr, seg): .expand_dims(seg=[seg]) ) - # ensure nans are only at the beginning - assert (this_surge_noadapt.notnull() | lslr_too_low).all(), seg + # ensure nans are only at the low end of LSLR or high end of rh_diff + rh_diff_too_high = ( + this_rh_diff_noadapt > this_surge_lookup.rh_diff.max() + ) + assert ( + this_surge_noadapt.notnull() | lslr_too_low | rh_diff_too_high + ).all(), seg surge_noadapt.append(this_surge_noadapt.fillna(0)) @@ -399,7 +416,7 @@ def _check_vals(rh_diff_arr, lslr_arr, seg): .reset_coords(drop=True) ) - # ensure nans are only at the beginning + # ensure nans are only at the low end of lslr assert (this_surge_adapt.notnull() | lslr_too_low).all(), seg surge_adapt.append(this_surge_adapt.fillna(0)) @@ -1500,9 +1517,11 @@ def get_refA( **params.refA_scenario_selectors, ) slr = slr.unstack("scen_mc") - slr = slr.squeeze(drop=True) + slr = slr.squeeze( + [d for d in slr.dims if len(slr[d]) == 1 and d != "seg"], drop=True + ) - costs, refA = calc_costs( + return calc_costs( inputs, slr, surge_lookup=surge, return_year0_hts=True, **model_kwargs ) From dc41e2d75a27cd301998a8ed6bf4df7460b13418 Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Tue, 7 Jan 2025 16:10:23 +0000 Subject: [PATCH 08/10] use DataArray.cumulative functionality in _get_lslr_plan_data --- pyCIAM/utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyCIAM/utils.py b/pyCIAM/utils.py index c96fee2..fd635b4 100644 --- a/pyCIAM/utils.py +++ b/pyCIAM/utils.py @@ -94,8 +94,7 @@ def _get_lslr_plan_data( if diaz_negative_retreat: RH_heights = _pos(RH_heights) else: - for i in range(1, len(RH_heights.at)): - RH_heights[{"at": i}] = RH_heights.isel(at=slice(None, i + 1)).max("at") + RH_heights = RH_heights.cumulative("at").max() # set initial RH_heights to 0 (e.g. # assuming no retreat or protection anywhere such that both w/ and w/o climate From 359d81515f561e0240c3cdfd084c81c00a6d196b Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Tue, 7 Jan 2025 16:30:24 +0000 Subject: [PATCH 09/10] undo temporary debugging return stmt --- pyCIAM/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index ef5018c..a42d18a 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -1521,7 +1521,7 @@ def get_refA( [d for d in slr.dims if len(slr[d]) == 1 and d != "seg"], drop=True ) - return calc_costs( + costs, refA = calc_costs( inputs, slr, surge_lookup=surge, return_year0_hts=True, **model_kwargs ) From 8f2505d427e96fd1005aede67c92b937af25c4de Mon Sep 17 00:00:00 2001 From: Ian Bolliger Date: Tue, 7 Jan 2025 16:39:58 +0000 Subject: [PATCH 10/10] bugfix in just-added surge lookup checks --- pyCIAM/run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyCIAM/run.py b/pyCIAM/run.py index a42d18a..873f911 100644 --- a/pyCIAM/run.py +++ b/pyCIAM/run.py @@ -401,7 +401,7 @@ def _check_vals(rh_diff_arr, lslr_arr, seg, check_rh_diff_max=True): surge_adapt = [] this_rh_diff_adapt = rh_diff.sel(seg=seg, drop=True) - _check_vals(this_rh_diff_noadapt, this_lslr, seg) + _check_vals(this_rh_diff_adapt, this_lslr, seg) for adapttype in this_surge_lookup.adapttype.values: this_surge_adapt = (