diff --git a/postprocessing.py b/postprocessing.py index 94dc3a7f..e9cb0bed 100644 --- a/postprocessing.py +++ b/postprocessing.py @@ -267,3 +267,48 @@ def merge_dicts(dict1, dict2): else: merged[key] = value return merged + + +# Stabilization constants for GCM change factors ("deltas") that are applied to +# historical baselines as ratios. The epsilon fraction sets the offset relative +# to the stream's historical mean flow; the floor keeps the offset nonzero for +# all-zero streams; the cap bounds how much any single statistic can be scaled. +RATIO_EPSILON_FRACTION = 0.01 +RATIO_EPSILON_FLOOR = 0.0001 +RATIO_CAP = 10.0 + + +def scale_aware_epsilon(values): + """ + Compute a stabilizing offset for ratio-based change factors, sized to the + magnitude of the historical data so that the same code works for large + rivers and small creeks alike. + Args: + values (list of float): Historical baseline values (e.g. all doy_mean + values for one model), used to establish the stream's flow scale + Returns: + float: RATIO_EPSILON_FRACTION times the mean of values, but never less + than RATIO_EPSILON_FLOOR + """ + if not values: + return RATIO_EPSILON_FLOOR + epsilon = RATIO_EPSILON_FRACTION * (sum(values) / len(values)) + return max(epsilon, RATIO_EPSILON_FLOOR) + + +def stabilized_ratio(projected, historical, epsilon): + """ + Ratio of projected to historical, stabilized for near-zero baselines. + Adding the same epsilon to numerator and denominator makes the ratio + approach 1 (no change) as both values approach zero, instead of exploding + when only the denominator is small. The result is clamped to + [1/RATIO_CAP, RATIO_CAP] as a backstop against noisy statistics. + Args: + projected (float): Future (projected) statistic value + historical (float): Historical statistic value for the same model + epsilon (float): Positive stabilizing offset from scale_aware_epsilon() + Returns: + float: Clamped change factor suitable for scaling a historical baseline + """ + ratio = (projected + epsilon) / (historical + epsilon) + return min(max(ratio, 1.0 / RATIO_CAP), RATIO_CAP) diff --git a/routes/arctic_hydrology.py b/routes/arctic_hydrology.py index ad9bd216..c86ec188 100644 --- a/routes/arctic_hydrology.py +++ b/routes/arctic_hydrology.py @@ -23,7 +23,11 @@ ) from fetch_data import fetch_data, fetch_layer_data, describe_via_wcps from validate_request import get_axis_encodings -from postprocessing import prune_nulls_with_max_intensity +from postprocessing import ( + prune_nulls_with_max_intensity, + scale_aware_epsilon, + stabilized_ratio, +) from csv_functions import create_csv from config import RAS_BASE_URL from . import routes @@ -490,6 +494,11 @@ def calculate_and_apply_gcm_diffs_to_blaskey_climatology(data_dict): """ Function to calculate GCM-projected changes in streamflow and apply them to the historical Blaskey climatology. Models without a '1990-2021' era (e.g. PGW models with no historical baseline) are silently skipped. + + Change factors are stabilized ratios (see postprocessing.stabilized_ratio): + a symmetric offset scaled to the model's historical mean flow keeps the + factor near 1 when both values are near zero (e.g. frozen winter minimums), + and factors are clamped to [1/RATIO_CAP, RATIO_CAP]. Args: data_dict (dict): Climatology data dict keyed by model then era Returns: @@ -502,6 +511,13 @@ def calculate_and_apply_gcm_diffs_to_blaskey_climatology(data_dict): continue if "1990-2021" not in data_dict[model]: continue + epsilon = scale_aware_epsilon( + [ + row["doy_mean"] + for row in data_dict[model]["1990-2021"] + if "doy_mean" in row + ] + ) adjusted_data_dict[model] = {} for era in data_dict[model].keys(): if era == "1990-2021": @@ -519,12 +535,20 @@ def calculate_and_apply_gcm_diffs_to_blaskey_climatology(data_dict): blaskey_historical = data_dict["historical"]["1990-2021"][i][stat] gcm_historical = data_dict[model]["1990-2021"][i][stat] gcm_projected = entry[stat] - denominator = gcm_historical - if denominator == 0: - denominator = 0.0001 - projected_quotient = gcm_projected / denominator + projected_quotient = stabilized_ratio( + gcm_projected, gcm_historical, epsilon + ) blaskey_adjusted = round(blaskey_historical * projected_quotient, 3) doy_stats[stat] = blaskey_adjusted + # each stat is scaled by its own ratio, so the adjusted values + # can cross; clamp to preserve min <= mean <= max + if all(k in doy_stats for k in ("doy_min", "doy_mean", "doy_max")): + doy_stats["doy_min"] = min( + doy_stats["doy_min"], doy_stats["doy_mean"] + ) + doy_stats["doy_max"] = max( + doy_stats["doy_max"], doy_stats["doy_mean"] + ) adjusted_data_dict[model][era].append(doy_stats) return adjusted_data_dict diff --git a/routes/conus_hydrology.py b/routes/conus_hydrology.py index 3631a5bf..c5f94ea2 100644 --- a/routes/conus_hydrology.py +++ b/routes/conus_hydrology.py @@ -25,7 +25,11 @@ ) from fetch_data import fetch_data, fetch_layer_data, describe_via_wcps from validate_request import get_axis_encodings -from postprocessing import prune_nulls_with_max_intensity +from postprocessing import ( + prune_nulls_with_max_intensity, + scale_aware_epsilon, + stabilized_ratio, +) from csv_functions import create_csv from config import RAS_BASE_URL from . import routes @@ -784,6 +788,11 @@ def calculate_and_apply_gcm_diffs_to_maurer_climatology(data_dict): Function to calculate the GCM-projected changes in streamflow stats and apply those changes to the historical Maurer climatology stats. This is done by first calculating the ratio between the GCM-projected future stat values and the GCM historical stat values, then applying that scaling factor to the Maurer historical stat values. + + Change factors are stabilized ratios (see postprocessing.stabilized_ratio): + a symmetric offset scaled to the model's historical mean flow keeps the + factor near 1 when both values are near zero (e.g. dry-season minimums), + and factors are clamped to [1/RATIO_CAP, RATIO_CAP]. Args: data_dict (dict): Data dictionary with the hydrology data populated Returns: @@ -795,6 +804,13 @@ def calculate_and_apply_gcm_diffs_to_maurer_climatology(data_dict): if model == "Maurer": adjusted_data_dict[model] = data_dict[model] continue + epsilon = scale_aware_epsilon( + [ + row["doy_mean"] + for row in data_dict[model]["historical"]["1976-2005"] + if "doy_mean" in row + ] + ) if model not in adjusted_data_dict: adjusted_data_dict[model] = {} for scenario in data_dict[model].keys(): @@ -819,14 +835,24 @@ def calculate_and_apply_gcm_diffs_to_maurer_climatology(data_dict): stat ] gcm_projected = entry[stat] - denominator = gcm_historical - if denominator == 0: - denominator = 0.0001 - projected_quotient = gcm_projected / denominator + projected_quotient = stabilized_ratio( + gcm_projected, gcm_historical, epsilon + ) maurer_adjusted = round( maurer_historical * projected_quotient, 3 ) doy_stats[stat] = maurer_adjusted + # each stat is scaled by its own ratio, so the adjusted + # values can cross; clamp to preserve min <= mean <= max + if all( + k in doy_stats for k in ("doy_min", "doy_mean", "doy_max") + ): + doy_stats["doy_min"] = min( + doy_stats["doy_min"], doy_stats["doy_mean"] + ) + doy_stats["doy_max"] = max( + doy_stats["doy_max"], doy_stats["doy_mean"] + ) adjusted_data_dict[model][scenario][era].append(doy_stats) return adjusted_data_dict @@ -962,6 +988,16 @@ def run_get_conus_hydrology_modeled_climatology(stream_id): data_dict = populate_feature_name_and_location_attributes(data_dict, gdf) data_dict = prune_nulls_with_max_intensity(data_dict) + # apply GCM-projected changes to Maurer climatology stats if source is "gcm_diff_applied_to_maurer" + # otherwise, if source is "original_gcm", then we are just returning the original GCM stats with no adjustments + if source == "gcm_diff_applied_to_maurer": + for landcover in data_dict["data"]: + data_dict["data"][landcover] = ( + calculate_and_apply_gcm_diffs_to_maurer_climatology( + data_dict["data"][landcover] + ) + ) + if request.args.get("format") == "csv": try: return create_csv( @@ -976,16 +1012,6 @@ def run_get_conus_hydrology_modeled_climatology(stream_id): except Exception as exc: return render_template("500/server_error.html"), 500 - # apply GCM-projected changes to Maurer climatology stats if source is "gcm_diff_applied_to_maurer" - # otherwise, if source is "original_gcm", then we are just returning the original GCM stats with no adjustments - if source == "gcm_diff_applied_to_maurer": - for landcover in data_dict["data"]: - data_dict["data"][landcover] = ( - calculate_and_apply_gcm_diffs_to_maurer_climatology( - data_dict["data"][landcover] - ) - ) - return jsonify(data_dict) except Exception as exc: diff --git a/tests/test_hydrology_adjustments.py b/tests/test_hydrology_adjustments.py new file mode 100644 index 00000000..188ef13c --- /dev/null +++ b/tests/test_hydrology_adjustments.py @@ -0,0 +1,233 @@ +import pytest +from postprocessing import ( + scale_aware_epsilon, + stabilized_ratio, + RATIO_CAP, + RATIO_EPSILON_FLOOR, + RATIO_EPSILON_FRACTION, +) +from routes.arctic_hydrology import ( + calculate_and_apply_gcm_diffs_to_blaskey_climatology, +) +from routes.conus_hydrology import ( + calculate_and_apply_gcm_diffs_to_maurer_climatology, +) + + +############################## +# 1. stabilized_ratio helper # +############################## + + +def test_stabilized_ratio_near_zero_baseline_does_not_explode(): + """ + Reproduces the production blowup at arctic stream 81009008, doy 119: + historical min 0.001 cfs vs projected min 492.363 cfs produced a raw + ratio of ~492,363x. With a scale-aware epsilon the factor must stay + within the cap. + """ + epsilon = scale_aware_epsilon([600.0] * 366) # ~600 cfs mean flow stream + ratio = stabilized_ratio(492.363, 0.001, epsilon) + assert ratio <= RATIO_CAP + + +def test_stabilized_ratio_zero_denominator_no_error(): + epsilon = scale_aware_epsilon([1.0]) + assert stabilized_ratio(5.0, 0.0, epsilon) <= RATIO_CAP + + +def test_stabilized_ratio_both_near_zero_approaches_one(): + epsilon = scale_aware_epsilon([100.0]) + ratio = stabilized_ratio(0.001, 0.002, epsilon) + assert 0.9 < ratio < 1.1 + + +def test_stabilized_ratio_ordinary_values_unchanged(): + """Healthy baselines should give (nearly) the plain ratio.""" + epsilon = scale_aware_epsilon([100.0]) + ratio = stabilized_ratio(200.0, 100.0, epsilon) + assert ratio == pytest.approx(2.0, rel=0.02) + + +def test_stabilized_ratio_clamps_low_end(): + epsilon = scale_aware_epsilon([100.0]) + assert stabilized_ratio(0.001, 5000.0, epsilon) == pytest.approx(1.0 / RATIO_CAP) + + +def test_scale_aware_epsilon_scales_with_flow(): + small = scale_aware_epsilon([1.0] * 10) + large = scale_aware_epsilon([10000.0] * 10) + assert small == pytest.approx(RATIO_EPSILON_FRACTION * 1.0) + assert large == pytest.approx(RATIO_EPSILON_FRACTION * 10000.0) + + +def test_scale_aware_epsilon_floor(): + assert scale_aware_epsilon([]) == RATIO_EPSILON_FLOOR + assert scale_aware_epsilon([0.0, 0.0]) == RATIO_EPSILON_FLOOR + + +########################################### +# 2. Arctic Blaskey climatology adjustment # +########################################### + + +def _arctic_data_dict(): + """ + Minimal one-day fixture mirroring production values at stream 81009008, + doy 119 (the observed 602,159.949 cfs doy_min blowup). + """ + day = {"doy": 119, "water_year_index": 211} + return { + "historical": { + "1990-2021": [ + dict(day, doy_min=1.223, doy_mean=3152.662, doy_max=10496.285) + ], + }, + "C2LE2": { + "1990-2021": [ + dict(day, doy_min=0.001, doy_mean=3964.098, doy_max=11532.644) + ], + "2034-2065": [ + dict(day, doy_min=492.363, doy_mean=3980.97, doy_max=9345.952) + ], + }, + # PGW model with no historical era: must be skipped, not crash + "PGWh": { + "2034-2065": [ + dict(day, doy_min=72.505, doy_mean=4257.796, doy_max=11770.358) + ], + }, + } + + +def test_blaskey_adjustment_min_stays_bounded(): + adjusted = calculate_and_apply_gcm_diffs_to_blaskey_climatology( + _arctic_data_dict() + ) + entry = adjusted["C2LE2"]["2034-2065"][0] + # blaskey min (1.223) can be scaled at most by RATIO_CAP + assert entry["doy_min"] <= 1.223 * RATIO_CAP + assert entry["doy_min"] < 602159.949 + + +def test_blaskey_adjustment_ordinary_stats_close_to_plain_ratio(): + adjusted = calculate_and_apply_gcm_diffs_to_blaskey_climatology( + _arctic_data_dict() + ) + entry = adjusted["C2LE2"]["2034-2065"][0] + plain_mean = 3152.662 * (3980.97 / 3964.098) + assert entry["doy_mean"] == pytest.approx(plain_mean, rel=0.05) + + +def test_blaskey_adjustment_skips_pgw_and_keeps_historical(): + adjusted = calculate_and_apply_gcm_diffs_to_blaskey_climatology( + _arctic_data_dict() + ) + assert "PGWh" not in adjusted + assert adjusted["historical"]["1990-2021"][0]["doy_min"] == 1.223 + assert "1990-2021" not in adjusted["C2LE2"] + + +########################################## +# 3. CONUS Maurer climatology adjustment # +########################################## + + +def _conus_data_dict(): + """One-day fixture in the CONUS shape (model -> scenario -> era).""" + day = {"doy": 119, "water_year_index": 211} + return { + "Maurer": { + "historical": { + "1976-2005": [ + dict(day, doy_min=296.9, doy_mean=708.793, doy_max=1373.0) + ], + }, + }, + "CCSM4": { + "historical": { + "1976-2005": [ + dict(day, doy_min=286.1, doy_mean=966.628, doy_max=4546.0) + ], + }, + "rcp45": { + "2016-2045": [ + dict(day, doy_min=0.002, doy_mean=900.0, doy_max=4000.0) + ], + }, + }, + } + + +def test_maurer_adjustment_near_zero_projected_min_stays_bounded(): + adjusted = calculate_and_apply_gcm_diffs_to_maurer_climatology(_conus_data_dict()) + entry = adjusted["CCSM4"]["rcp45"]["2016-2045"][0] + # projected min ~0 against healthy historical min: clamped at 1/RATIO_CAP + assert entry["doy_min"] >= 296.9 / RATIO_CAP - 0.001 + assert entry["doy_min"] <= 296.9 + + +def test_maurer_adjustment_preserves_structure(): + adjusted = calculate_and_apply_gcm_diffs_to_maurer_climatology(_conus_data_dict()) + assert adjusted["Maurer"]["historical"]["1976-2005"][0]["doy_mean"] == 708.793 + assert "rcp45" in adjusted["CCSM4"] + + +#################################### +# 4. min <= mean <= max invariant # +#################################### + + +def _assert_ordered(entry): + assert entry["doy_min"] <= entry["doy_mean"] <= entry["doy_max"] + + +def test_blaskey_adjustment_enforces_stat_ordering(): + """ + Each stat is scaled by its own ratio, so a large projected minimum can + overtake the adjusted mean. The clamp must restore ordering. + """ + day = {"doy": 200, "water_year_index": 292} + data = { + "historical": { + "1990-2021": [dict(day, doy_min=100.0, doy_mean=200.0, doy_max=1000.0)], + }, + "MODEL": { + "1990-2021": [dict(day, doy_min=100.0, doy_mean=200.0, doy_max=4000.0)], + "2034-2065": [dict(day, doy_min=400.0, doy_mean=200.0, doy_max=1200.0)], + }, + } + adjusted = calculate_and_apply_gcm_diffs_to_blaskey_climatology(data) + _assert_ordered(adjusted["MODEL"]["2034-2065"][0]) + + +def test_maurer_adjustment_enforces_stat_ordering(): + """ + Production values from CONUS stream 50563 (dynamic / CCSM4 / rcp26, + doy 214), which yielded adjusted doy_min 353.042 > doy_max 343.703 + because the GCM historical max held a flood day that does not recur. + """ + day = {"doy": 214, "water_year_index": 306} + data = { + "Maurer": { + "historical": { + "1976-2005": [ + dict(day, doy_min=296.9, doy_mean=708.793, doy_max=1373.0) + ], + }, + }, + "CCSM4": { + "historical": { + "1976-2005": [ + dict(day, doy_min=286.1, doy_mean=966.628, doy_max=4546.0) + ], + }, + "rcp26": { + "2016-2045": [ + dict(day, doy_min=340.2, doy_mean=627.062, doy_max=1138.0) + ], + }, + }, + } + adjusted = calculate_and_apply_gcm_diffs_to_maurer_climatology(data) + _assert_ordered(adjusted["CCSM4"]["rcp26"]["2016-2045"][0])