From c0932ba66d66b565f912ff2ec3f105250c7bc248 Mon Sep 17 00:00:00 2001 From: btobers Date: Thu, 14 Aug 2025 16:05:08 -0400 Subject: [PATCH 1/4] return precip bias adjustment factors to run_simulation, warn on factor greater than 2x --- pygem/bin/run/run_simulation.py | 41 ++++++++++++++++++++++----------- pygem/gcmbiasadj.py | 14 ++++------- 2 files changed, 31 insertions(+), 24 deletions(-) diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index df3ba04a..8237cc94 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -23,6 +23,7 @@ import os import sys import time +import warnings import matplotlib.pyplot as plt import numpy as np @@ -457,6 +458,7 @@ def run(list_packed_vars): gcm_elev_adj = gcm_elev gcm_temp_adj = gcm_temp[:, sim_idx_start:] gcm_prec_adj = gcm_prec[:, sim_idx_start:] + gcm_prec_biasadj_frac = np.ones(gcm_prec_adj.shape[0]) # Bias correct based on reference climate data else: # OPTION 1: Adjust temp using Huss and Hock (2015), prec similar but addresses for variance and outliers @@ -472,14 +474,16 @@ def run(list_packed_vars): args.ref_startyear, ) # Precipitation bias correction - gcm_prec_adj, gcm_elev_adj = gcmbiasadj.prec_biasadj_opt1( - ref_prec, - ref_elev, - gcm_prec, - dates_table_ref, - dates_table_full, - args.sim_startyear, - args.ref_startyear, + gcm_prec_adj, gcm_elev_adj, gcm_prec_biasadj_frac = ( + gcmbiasadj.prec_biasadj_opt1( + ref_prec, + ref_elev, + gcm_prec, + dates_table_ref, + dates_table_full, + args.sim_startyear, + args.ref_startyear, + ) ) # OPTION 2: Adjust temp and prec using Huss and Hock (2015) elif args.option_bias_adjustment == 2: @@ -494,12 +498,14 @@ def run(list_packed_vars): args.ref_startyear, ) # Precipitation bias correction - gcm_prec_adj, gcm_elev_adj = gcmbiasadj.prec_biasadj_HH2015( - ref_prec, - ref_elev, - gcm_prec, - dates_table_ref, - dates_table_full, + gcm_prec_adj, gcm_elev_adj, gcm_prec_biasadj_frac = ( + gcmbiasadj.prec_biasadj_HH2015( + ref_prec, + ref_elev, + gcm_prec, + dates_table_ref, + dates_table_full, + ) ) # OPTION 3: Adjust temp and prec using quantile delta mapping, Cannon et al. (2015) elif args.option_bias_adjustment == 3: @@ -524,6 +530,7 @@ def run(list_packed_vars): args.sim_startyear, args.ref_startyear, ) + gcm_prec_biasadj_frac = np.ones(gcm_prec_adj.shape[0]) # assert that the gcm_elev_adj is not None assert gcm_elev_adj is not None, 'No GCM elevation data' @@ -640,6 +647,12 @@ def run(list_packed_vars): 'prec': gcm_prec_adj[glac, :], 'lr': gcm_lr[glac, :], } + fact = gcm_prec_biasadj_frac[glac] + if fact < 0.5 or fact > 2: + warnings.warn( + f'GCM precipitation bias adustment factor for {glacier_str} greater than 2 ({round(fact, 2)})', + Warning, + ) gdir.dates_table = dates_table glacier_area_km2 = fls[0].widths_m * fls[0].dx_meter / 1e6 diff --git a/pygem/gcmbiasadj.py b/pygem/gcmbiasadj.py index c894ea2e..eaaa3db6 100755 --- a/pygem/gcmbiasadj.py +++ b/pygem/gcmbiasadj.py @@ -298,12 +298,9 @@ def prec_biasadj_HH2015( gcm_prec_biasadj_subset = gcm_prec_biasadj[ :, gcm_subset_idx_start : gcm_subset_idx_end + 1 ][:, gcm_spinupyears * 12 :] - gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / ref_prec_nospinup.sum( + gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / gcm_prec_nospinup.sum( axis=1 ) - assert np.min(gcm_prec_biasadj_frac) > 0.5 and np.max(gcm_prec_biasadj_frac) < 2, ( - 'Error with gcm precipitation bias adjustment: total ref and gcm prec differ by more than factor of 2' - ) assert gcm_prec_biasadj.max() <= 10, ( 'gcm_prec_adj (precipitation bias adjustment) too high, needs to be modified' ) @@ -311,7 +308,7 @@ def prec_biasadj_HH2015( 'gcm_prec_adj is producing a negative precipitation value' ) - return gcm_prec_biasadj, gcm_elev_biasadj + return gcm_prec_biasadj, gcm_elev_biasadj, gcm_prec_biasadj_frac def prec_biasadj_opt1( @@ -454,12 +451,9 @@ def prec_biasadj_opt1( gcm_prec_biasadj_subset = gcm_prec_biasadj[ :, gcm_subset_idx_start : gcm_subset_idx_end + 1 ][:, gcm_spinupyears * 12 :] - gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / ref_prec_nospinup.sum( + gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / gcm_prec_nospinup.sum( axis=1 ) - assert np.min(gcm_prec_biasadj_frac) > 0.5 and np.max(gcm_prec_biasadj_frac) < 2, ( - 'Error with gcm precipitation bias adjustment: total ref and gcm prec differ by more than factor of 2' - ) assert gcm_prec_biasadj.max() <= 10, ( 'gcm_prec_adj (precipitation bias adjustment) too high, needs to be modified' ) @@ -467,7 +461,7 @@ def prec_biasadj_opt1( 'gcm_prec_adj is producing a negative precipitation value' ) - return gcm_prec_biasadj, gcm_elev_biasadj + return gcm_prec_biasadj, gcm_elev_biasadj, gcm_prec_biasadj_frac def temp_biasadj_QDM( From 0458f726afb37d46b1a504c271abe204ac32e78e Mon Sep 17 00:00:00 2001 From: btobers Date: Thu, 14 Aug 2025 16:20:11 -0400 Subject: [PATCH 2/4] stacklevel for warnings.warn() --- pygem/bin/run/run_simulation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index 8237cc94..918ab3fa 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -652,6 +652,7 @@ def run(list_packed_vars): warnings.warn( f'GCM precipitation bias adustment factor for {glacier_str} greater than 2 ({round(fact, 2)})', Warning, + stacklevel=2, ) gdir.dates_table = dates_table From a4769b52962d8a42b521576e0b431a20bb51b31f Mon Sep 17 00:00:00 2001 From: btobers Date: Mon, 18 Aug 2025 10:05:04 -0400 Subject: [PATCH 3/4] bug fix in import of pygem config --- pygem/bin/op/compress_gdirs.py | 8 ++++---- pygem/bin/run/run_inversion.py | 8 ++++---- pygem/bin/run/run_spinup.py | 8 ++++---- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/pygem/bin/op/compress_gdirs.py b/pygem/bin/op/compress_gdirs.py index 9ada5cdb..7d121817 100644 --- a/pygem/bin/op/compress_gdirs.py +++ b/pygem/bin/op/compress_gdirs.py @@ -16,12 +16,12 @@ from oggm import utils, workflow # pygem imports -import pygem.setup.config as config +from pygem.setup.config import ConfigManager -# check for config -config.ensure_config() +# instantiate ConfigManager +config_manager = ConfigManager() # read the config -pygem_prms = config.read_config() +pygem_prms = config_manager.read_config() # Initialize OGGM subprocess cfg.initialize(logging_level='WARNING') diff --git a/pygem/bin/run/run_inversion.py b/pygem/bin/run/run_inversion.py index b022f160..d118ad5a 100644 --- a/pygem/bin/run/run_inversion.py +++ b/pygem/bin/run/run_inversion.py @@ -6,12 +6,12 @@ import pandas as pd # pygem imports -import pygem.setup.config as config +from pygem.setup.config import ConfigManager -# check for config -config.ensure_config() +# instantiate ConfigManager +config_manager = ConfigManager() # read the config -pygem_prms = config.read_config() +pygem_prms = config_manager.read_config() from oggm import cfg, tasks, workflow import pygem.pygem_modelsetup as modelsetup diff --git a/pygem/bin/run/run_spinup.py b/pygem/bin/run/run_spinup.py index 182323bd..2ae7eb24 100644 --- a/pygem/bin/run/run_spinup.py +++ b/pygem/bin/run/run_spinup.py @@ -7,12 +7,12 @@ import pandas as pd # pygem imports -import pygem.setup.config as config +from pygem.setup.config import ConfigManager -# check for config -config.ensure_config() +# instantiate ConfigManager +config_manager = ConfigManager() # read the config -pygem_prms = config.read_config() +pygem_prms = config_manager.read_config() from oggm import tasks, workflow import pygem.pygem_modelsetup as modelsetup From 23a0a9d9638faaedce3b30dab9be8f07f4b606eb Mon Sep 17 00:00:00 2001 From: btobers Date: Mon, 18 Aug 2025 12:42:48 -0400 Subject: [PATCH 4/4] remove spinup years, revert gcm_prec_biasadj_frac to gcm_prec_biasadj_subset/ref_prec --- pygem/bin/run/run_simulation.py | 2 +- pygem/gcmbiasadj.py | 99 +++++++++------------------------ 2 files changed, 27 insertions(+), 74 deletions(-) diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index 918ab3fa..1225bd50 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -650,7 +650,7 @@ def run(list_packed_vars): fact = gcm_prec_biasadj_frac[glac] if fact < 0.5 or fact > 2: warnings.warn( - f'GCM precipitation bias adustment factor for {glacier_str} greater than 2 ({round(fact, 2)})', + f'Bias-adjusted GCM precipitation for {glacier_str} differs from that of the refernce climate data by a factor greater than 2x ({round(fact, 2)})', Warning, stacklevel=2, ) diff --git a/pygem/gcmbiasadj.py b/pygem/gcmbiasadj.py index eaaa3db6..89086dfa 100755 --- a/pygem/gcmbiasadj.py +++ b/pygem/gcmbiasadj.py @@ -75,8 +75,6 @@ def temp_biasadj_HH2015( dates_table, sim_startyear, ref_startyear, - ref_spinupyears=0, - gcm_spinupyears=0, debug=False, ): """ @@ -113,28 +111,20 @@ def temp_biasadj_HH2015( )[0][0] gcm_temp_subset = gcm_temp[:, gcm_subset_idx_start : gcm_subset_idx_end + 1] - # Remove spinup years, so adjustment performed over calibration period - ref_temp_nospinup = ref_temp[:, ref_spinupyears * 12 :] - gcm_temp_nospinup = gcm_temp_subset[:, gcm_spinupyears * 12 :] - # Roll months so they are aligned with simulation months roll_amt = -1 * (12 - gcm_subset_idx_start % 12) if roll_amt == -12: roll_amt = 0 # Mean monthly temperature - ref_temp_monthly_avg = np.roll( - monthly_avg_2darray(ref_temp_nospinup), roll_amt, axis=1 - ) + ref_temp_monthly_avg = np.roll(monthly_avg_2darray(ref_temp), roll_amt, axis=1) gcm_temp_monthly_avg = np.roll( - monthly_avg_2darray(gcm_temp_nospinup), roll_amt, axis=1 + monthly_avg_2darray(gcm_temp_subset), roll_amt, axis=1 ) # Standard deviation monthly temperature - ref_temp_monthly_std = np.roll( - monthly_std_2darray(ref_temp_nospinup), roll_amt, axis=1 - ) + ref_temp_monthly_std = np.roll(monthly_std_2darray(ref_temp), roll_amt, axis=1) gcm_temp_monthly_std = np.roll( - monthly_std_2darray(gcm_temp_nospinup), roll_amt, axis=1 + monthly_std_2darray(gcm_temp_subset), roll_amt, axis=1 ) # Monthly bias adjustment (additive) @@ -182,21 +172,17 @@ def temp_biasadj_HH2015( # Assert that mean temperatures for all the glaciers must be more-or-less equal gcm_temp_biasadj_subset = gcm_temp_biasadj[ :, gcm_subset_idx_start : gcm_subset_idx_end + 1 - ][:, ref_spinupyears * 12 :] + ] if sim_startyear == ref_startyear: if debug: print( - ( - np.mean(gcm_temp_biasadj_subset, axis=1) - - np.mean(ref_temp[:, ref_spinupyears * 12 :], axis=1) - ) + (np.mean(gcm_temp_biasadj_subset, axis=1) - np.mean(ref_temp, axis=1)) ) assert ( np.max( np.abs( - np.mean(gcm_temp_biasadj_subset, axis=1) - - np.mean(ref_temp[:, ref_spinupyears * 12 :], axis=1) + np.mean(gcm_temp_biasadj_subset, axis=1) - np.mean(ref_temp, axis=1) ) ) < 1 @@ -206,10 +192,7 @@ def temp_biasadj_HH2015( else: if debug: print( - ( - np.mean(gcm_temp_biasadj_subset, axis=1) - - np.mean(ref_temp[:, ref_spinupyears * 12 :], axis=1) - ) + (np.mean(gcm_temp_biasadj_subset, axis=1) - np.mean(ref_temp, axis=1)) ) return gcm_temp_biasadj, gcm_elev_biasadj @@ -223,8 +206,6 @@ def prec_biasadj_HH2015( dates_table, sim_startyear, ref_startyear, - ref_spinupyears=0, - gcm_spinupyears=0, ): """ Huss and Hock (2015) precipitation bias correction based on mean (multiplicative) @@ -256,20 +237,14 @@ def prec_biasadj_HH2015( )[0][0] gcm_prec_subset = gcm_prec[:, gcm_subset_idx_start : gcm_subset_idx_end + 1] - # Remove spinup years, so adjustment performed over calibration period - ref_prec_nospinup = ref_prec[:, ref_spinupyears * 12 :] - gcm_prec_nospinup = gcm_prec_subset[:, gcm_spinupyears * 12 :] - # Roll months so they are aligned with simulation months roll_amt = -1 * (12 - gcm_subset_idx_start % 12) # PRECIPITATION BIAS CORRECTIONS # Monthly mean precipitation - ref_prec_monthly_avg = np.roll( - monthly_avg_2darray(ref_prec_nospinup), roll_amt, axis=1 - ) + ref_prec_monthly_avg = np.roll(monthly_avg_2darray(ref_prec), roll_amt, axis=1) gcm_prec_monthly_avg = np.roll( - monthly_avg_2darray(gcm_prec_nospinup), roll_amt, axis=1 + monthly_avg_2darray(gcm_prec_subset), roll_amt, axis=1 ) bias_adj_prec_monthly = ref_prec_monthly_avg / gcm_prec_monthly_avg @@ -297,10 +272,8 @@ def prec_biasadj_HH2015( # Assertion that bias adjustment does not drastically modify the precipitation and are reasonable gcm_prec_biasadj_subset = gcm_prec_biasadj[ :, gcm_subset_idx_start : gcm_subset_idx_end + 1 - ][:, gcm_spinupyears * 12 :] - gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / gcm_prec_nospinup.sum( - axis=1 - ) + ] + gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / ref_prec.sum(axis=1) assert gcm_prec_biasadj.max() <= 10, ( 'gcm_prec_adj (precipitation bias adjustment) too high, needs to be modified' ) @@ -319,8 +292,6 @@ def prec_biasadj_opt1( dates_table, sim_startyear, ref_startyear, - ref_spinupyears=0, - gcm_spinupyears=0, ): """ Precipitation bias correction based on mean with limited maximum @@ -352,20 +323,14 @@ def prec_biasadj_opt1( )[0][0] gcm_prec_subset = gcm_prec[:, gcm_subset_idx_start : gcm_subset_idx_end + 1] - # Remove spinup years, so adjustment performed over calibration period - ref_prec_nospinup = ref_prec[:, ref_spinupyears * 12 :] - gcm_prec_nospinup = gcm_prec_subset[:, gcm_spinupyears * 12 :] - # Roll months so they are aligned with simulation months roll_amt = -1 * (12 - gcm_subset_idx_start % 12) # PRECIPITATION BIAS CORRECTIONS # Monthly mean precipitation - ref_prec_monthly_avg = np.roll( - monthly_avg_2darray(ref_prec_nospinup), roll_amt, axis=1 - ) + ref_prec_monthly_avg = np.roll(monthly_avg_2darray(ref_prec), roll_amt, axis=1) gcm_prec_monthly_avg = np.roll( - monthly_avg_2darray(gcm_prec_nospinup), roll_amt, axis=1 + monthly_avg_2darray(gcm_prec_subset), roll_amt, axis=1 ) bias_adj_prec_monthly = ref_prec_monthly_avg / gcm_prec_monthly_avg @@ -388,9 +353,7 @@ def prec_biasadj_opt1( ) # Adjust variance based on zscore and reference standard deviation - ref_prec_monthly_std = np.roll( - monthly_std_2darray(ref_prec_nospinup), roll_amt, axis=1 - ) + ref_prec_monthly_std = np.roll(monthly_std_2darray(ref_prec), roll_amt, axis=1) gcm_prec_biasadj_raw_monthly_avg = monthly_avg_2darray( gcm_prec_biasadj_raw[:, 0 : ref_prec.shape[1]] ) @@ -412,9 +375,9 @@ def prec_biasadj_opt1( # Identify outliers using reference's monthly maximum adjusted for future increases ref_prec_monthly_max = np.roll( ( - ref_prec_nospinup.reshape(-1, 12) + ref_prec.reshape(-1, 12) .transpose() - .reshape(-1, int(ref_prec_nospinup.shape[1] / 12)) + .reshape(-1, int(ref_prec.shape[1] / 12)) .max(1) .reshape(12, -1) .transpose() @@ -450,10 +413,8 @@ def prec_biasadj_opt1( # Assertion that bias adjustment does not drastically modify the precipitation and are reasonable gcm_prec_biasadj_subset = gcm_prec_biasadj[ :, gcm_subset_idx_start : gcm_subset_idx_end + 1 - ][:, gcm_spinupyears * 12 :] - gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / gcm_prec_nospinup.sum( - axis=1 - ) + ] + gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / ref_prec.sum(axis=1) assert gcm_prec_biasadj.max() <= 10, ( 'gcm_prec_adj (precipitation bias adjustment) too high, needs to be modified' ) @@ -472,8 +433,6 @@ def temp_biasadj_QDM( dates_table, sim_startyear, ref_startyear, - ref_spinupyears=0, - gcm_spinupyears=0, ): """ Cannon et al. (2015) temperature bias correction based on quantile delta mapping @@ -518,9 +477,9 @@ def temp_biasadj_QDM( )[0][0] gcm_temp_historic = gcm_temp[:, gcm_subset_idx_start : gcm_subset_idx_end + 1] - # Remove spinup years, so adjustment performed over calibration period - ref_temp_nospinup = ref_temp[:, ref_spinupyears * 12 :] + 273.15 - gcm_temp_nospinup = gcm_temp_historic[:, gcm_spinupyears * 12 :] + 273.15 + # Convert to Kelvin + ref_temp = ref_temp + 273.15 + gcm_temp_historic = gcm_temp_historic + 273.15 # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), @@ -564,8 +523,8 @@ def temp_biasadj_QDM( for ival, projected_value in enumerate(bc_temp_loop): percentile = percentileofscore(bc_temp_loop, projected_value) bias_correction_factor = np.percentile( - ref_temp_nospinup, percentile - ) / np.percentile(gcm_temp_nospinup, percentile) + ref_temp, percentile + ) / np.percentile(gcm_temp_historic, percentile) bc_temp_loop_corrected[ival] = projected_value * bias_correction_factor # append the values from each time period to a list gcm_temp_biasadj.append(bc_temp_loop_corrected) @@ -596,8 +555,6 @@ def prec_biasadj_QDM( dates_table, sim_startyear, ref_startyear, - ref_spinupyears=0, - gcm_spinupyears=0, ): """ Cannon et al. (2015) precipitation bias correction based on quantile delta mapping @@ -643,10 +600,6 @@ def prec_biasadj_QDM( )[0][0] gcm_prec_historic = gcm_prec[:, gcm_subset_idx_start : gcm_subset_idx_end + 1] - # Remove spinup years, so adjustment performed over calibration period - ref_prec_nospinup = ref_prec[:, ref_spinupyears * 12 :] - gcm_prec_nospinup = gcm_prec_historic[:, gcm_spinupyears * 12 :] - # if/else statement for whether or not the full GCM period is the same as the simulation period # create GCM subset for applying bias-correction (e.g., 2000-2100), # that does not include the earlier reference years (e.g., 1981-2000) @@ -687,8 +640,8 @@ def prec_biasadj_QDM( for ival, projected_value in enumerate(bc_prec_loop): percentile = percentileofscore(bc_prec_loop, projected_value) bias_correction_factor = np.percentile( - ref_prec_nospinup, percentile - ) / np.percentile(gcm_prec_nospinup, percentile) + ref_prec, percentile + ) / np.percentile(gcm_prec_historic, percentile) bc_prec_loop_corrected[ival] = projected_value * bias_correction_factor # append the values from each time period to a list gcm_prec_biasadj.append(bc_prec_loop_corrected)