From a82fb97b90f38f2b465686508cf0705215714c3e Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 26 Jun 2023 09:48:35 +0200 Subject: [PATCH 1/3] added more flexiblity to dynamic spinup --- oggm/core/dynamic_spinup.py | 197 ++++++++++++++++++++---------------- oggm/tests/test_models.py | 57 +++++++++-- 2 files changed, 155 insertions(+), 99 deletions(-) diff --git a/oggm/core/dynamic_spinup.py b/oggm/core/dynamic_spinup.py index 754fbdac0..a2f2bb632 100644 --- a/oggm/core/dynamic_spinup.py +++ b/oggm/core/dynamic_spinup.py @@ -34,7 +34,8 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, mb_model_historical=None, mb_model_spinup=None, spinup_period=20, spinup_start_yr=None, min_spinup_period=10, spinup_start_yr_max=None, - yr_rgi=None, minimise_for='area', precision_percent=1, + yr_target=None, target_value=None, + minimise_for='area', precision_percent=1, precision_absolute=1, min_ice_thickness=None, first_guess_t_bias=-2, t_bias_max_step_length=2, maxiter=30, output_filesuffix='_dynamic_spinup', @@ -81,10 +82,10 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, period of spinup_start_yr until rgi_year (e.g. 1979 - 2000). spinup_period : int The period how long the spinup should run. Start date of historical run - is defined "yr_rgi - spinup_period". Minimum allowed value is 10. If + is defined "yr_target - spinup_period". Minimum allowed value is 10. If the provided climate data starts at year later than - (yr_rgi - spinup_period) the spinup_period is set to - (yr_rgi - yr_climate_start). Caution if spinup_start_yr is set the + (yr_target - spinup_period) the spinup_period is set to + (yr_target - yr_climate_start). Caution if spinup_start_yr is set the spinup_period is ignored. Default is 20 spinup_start_yr : int or None @@ -100,13 +101,20 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, spinup_start_yr_max : int or None Possibility to provide a maximum year where the dynamic spinup must start from at least. If set, this overrides the min_spinup_period if - yr_rgi - spinup_start_yr_max > min_spinup_period. + yr_target - spinup_start_yr_max > min_spinup_period. Default is None - yr_rgi : int - The rgi date, at which we want to match area or volume. + yr_target : int or None + The year at which we want to match area or volume. If None, gdir.rgi_date + 1 is used (the default). + Default is None + target_value : float or None + The value we want to match at yr_target. Depending on minimise_for this + value is interpreted as an area in km2 or a volume in km3. If None the + total area or volume from the provided initial flowlines is used. + Default is None minimise_for : str - The variable we want to match at yr_rgi. Options are 'area' or 'volume'. + The variable we want to match at yr_target. Options are 'area' or + 'volume'. Default is 'area' precision_percent : float Gives the precision we want to match in percent. The algorithm makes @@ -135,8 +143,8 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, MassBalanceModel in °C. Default is -2. t_bias_max_step_length : float - Defines the maximums allowed change of t_bias between two iterations. Is - needed to avoid to large changes. + Defines the maximums allowed change of t_bias between two iterations. + Is needed to avoid to large changes. Default is 2. maxiter : int Maximum number of minimisation iterations per spinup period. If reached @@ -169,11 +177,11 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, the returned value is np.nan. Default is False ye : int - end year of the model run, must be larger than yr_rgi. If nothing is - given it is set to yr_rgi. It is not recommended to use it if only data - until yr_rgi is needed for calibration as this increases the run time - of each iteration during the iterative minimisation. Instead use - run_from_climate_data afterwards and merge both outputs using + end year of the model run, must be larger than yr_target. If nothing is + given it is set to yr_target. It is not recommended to use it if only + data until yr_target is needed for calibration as this increases the + run time of each iteration during the iterative minimisation. Instead + use run_from_climate_data afterwards and merge both outputs using merge_consecutive_run_outputs. Default is None model_flowline_filesuffix : str @@ -199,17 +207,17 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, evolution_model = decide_evolution_model(evolution_model) - if yr_rgi is None: + if yr_target is None: # Even in calendar dates, we prefer to set rgi_year in the next year # as the rgi is often from snow free images the year before (e.g. Aug) - yr_rgi = gdir.rgi_date + 1 + yr_target = gdir.rgi_date + 1 if ye is None: - ye = yr_rgi + ye = yr_target - if ye < yr_rgi: + if ye < yr_target: raise RuntimeError(f'The provided end year (ye = {ye}) must be larger' - f'than the rgi date (yr_rgi = {yr_rgi}!') + f'than the target year (yr_target = {yr_target}!') yr_min = gdir.get_climate_info()['baseline_yr_0'] @@ -229,8 +237,8 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, f'{spinup_start_yr}) must be smaller than ' f'the maximum start year ' f'{spinup_start_yr_max}!') - if (yr_rgi - spinup_start_yr_max) > min_spinup_period: - min_spinup_period = (yr_rgi - spinup_start_yr_max) + if (yr_target - spinup_start_yr_max) > min_spinup_period: + min_spinup_period = (yr_target - spinup_start_yr_max) if init_model_filesuffix is not None: fp = gdir.get_filepath('model_geometry', @@ -247,7 +255,7 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, else: fls_spinup = copy.deepcopy(init_model_fls) - # MassBalance for actual run from yr_spinup to yr_rgi + # MassBalance for actual run from yr_spinup to yr_target if mb_model_historical is None: mb_model_historical = MultipleFlowlineMassBalance( gdir, mb_model_class=MonthlyTIModel, @@ -298,10 +306,10 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None, # this function saves a model without conducting a dynamic spinup, but with # the provided output_filesuffix, so following tasks can find it. - # This is necessary if yr_rgi < yr_min + 10 or if the dynamic spinup failed. + # This is necessary if yr_target < yr_min + 10 or if the dynamic spinup failed. def save_model_without_dynamic_spinup(): gdir.add_to_diagnostics('run_dynamic_spinup_success', False) - yr_use = np.clip(yr_rgi, yr_min, None) + yr_use = np.clip(yr_target, yr_min, None) model_dynamic_spinup_end = evolution_model(fls_spinup, mb_model_historical, y0=yr_use, @@ -322,12 +330,12 @@ def save_model_without_dynamic_spinup(): return model_dynamic_spinup_end - if yr_rgi < yr_min + min_spinup_period: + if yr_target < yr_min + min_spinup_period: log.warning('The provided rgi_date is smaller than yr_climate_start + ' 'min_spinup_period, therefore no dynamic spinup is ' 'conducted and the original flowlines are saved at the ' - 'provided rgi_date or the start year of the provided ' - 'climate data (if yr_climate_start > yr_rgi)') + 'provided target year or the start year of the provided ' + 'climate data (if yr_climate_start > yr_target)') if ignore_errors: model_dynamic_spinup_end = save_model_without_dynamic_spinup() if return_t_bias_best: @@ -355,7 +363,10 @@ def save_model_without_dynamic_spinup(): else: raise NotImplementedError cost_var = f'{minimise_for}_{unit}' - reference_value = np.sum([getattr(f, cost_var) for f in fls_ref]) + if target_value is None: + reference_value = np.sum([getattr(f, cost_var) for f in fls_ref]) + else: + reference_value = target_value other_reference_value = np.sum([getattr(f, f'{other_variable}_{other_unit}') for f in fls_ref]) @@ -380,7 +391,7 @@ def save_model_without_dynamic_spinup(): forward_model_runs = [0] # the actual spinup run - def run_model_with_spinup_to_rgi_date(t_bias): + def run_model_with_spinup_to_target_year(t_bias): forward_model_runs.append(forward_model_runs[-1] + 1) # with t_bias the glacier state after spinup is changed between iterations @@ -425,11 +436,11 @@ def run_model_with_spinup_to_rgi_date(t_bias): if type(ds) == tuple: ds = ds[0] - model_area_km2 = ds.area_m2_min_h.loc[yr_rgi].values * 1e-6 - model_volume_km3 = ds.volume_m3.loc[yr_rgi].values * 1e-9 + model_area_km2 = ds.area_m2_min_h.loc[yr_target].values * 1e-6 + model_volume_km3 = ds.volume_m3.loc[yr_target].values * 1e-9 else: # only run to rgi date and extract values - model_historical.run_until(yr_rgi) + model_historical.run_until(yr_target) fls = model_historical.fls model_area_km2 = np.sum( [np.sum(fl.bin_area_m2[fl.thick > min_ice_thickness]) @@ -448,7 +459,7 @@ def run_model_with_spinup_to_rgi_date(t_bias): def cost_fct(t_bias, model_dynamic_spinup_end_loc, other_variable_mismatch_loc): # actual model run model_value, other_value, model_dynamic_spinup, ice_free = \ - run_model_with_spinup_to_rgi_date(t_bias) + run_model_with_spinup_to_target_year(t_bias) # save the final model for later model_dynamic_spinup_end_loc.append(copy.deepcopy(model_dynamic_spinup)) @@ -773,20 +784,22 @@ def get_mismatch(t_bias): # years and the second try is to use a period of 'min_spinup_period' years, # if it still fails the actual error is raised) if spinup_start_yr is not None: - spinup_period_initial = min(yr_rgi - spinup_start_yr, yr_rgi - yr_min) + spinup_period_initial = min(yr_target - spinup_start_yr, + yr_target - yr_min) else: - spinup_period_initial = min(spinup_period, yr_rgi - yr_min) + spinup_period_initial = min(spinup_period, yr_target - yr_min) if spinup_period_initial <= min_spinup_period: spinup_periods_to_try = [min_spinup_period] else: # try out a maximum of three different spinup_periods spinup_periods_to_try = [spinup_period_initial, - int((spinup_period_initial + min_spinup_period) / 2), + int((spinup_period_initial + + min_spinup_period) / 2), min_spinup_period] # after defining the initial spinup period we can define the year for the # fixed_geometry_spinup if add_fixed_geometry_spinup: - fixed_geometry_spinup_yr = yr_rgi - spinup_period_initial + fixed_geometry_spinup_yr = yr_target - spinup_period_initial else: fixed_geometry_spinup_yr = None @@ -797,14 +810,14 @@ def get_mismatch(t_bias): provided_mb_model_spinup = True for spinup_period in spinup_periods_to_try: - yr_spinup = yr_rgi - spinup_period + yr_spinup = yr_target - spinup_period if not provided_mb_model_spinup: # define spinup MassBalance - # spinup is running for 'yr_rgi - yr_spinup' years, using a + # spinup is running for 'yr_target - yr_spinup' years, using a # ConstantMassBalance - y0_spinup = (yr_spinup + yr_rgi) / 2 - halfsize_spinup = yr_rgi - y0_spinup + y0_spinup = (yr_spinup + yr_target) / 2 + halfsize_spinup = yr_target - y0_spinup mb_model_spinup = MultipleFlowlineMassBalance( gdir, mb_model_class=ConstantMassBalance, filename='climate_historical', @@ -871,7 +884,7 @@ def get_mismatch(t_bias): # For operational runs we ignore the warnings warnings.filterwarnings('ignore', category=RuntimeWarning) model_dynamic_spinup_end[-1].run_until_and_store( - yr_rgi, + yr_target, geom_path=geom_path, diag_path=diag_path, fl_diag_path=fl_diag_path, ) @@ -914,7 +927,7 @@ def dynamic_melt_f_run_with_dynamic_spinup( output_filesuffix='', evolution_model=None, mb_model_historical=None, mb_model_spinup=None, minimise_for='area', climate_input_filesuffix='', spinup_period=20, - min_spinup_period=10, yr_rgi=None, precision_percent=1, + min_spinup_period=10, yr_target=None, precision_percent=1, precision_absolute=1, min_ice_thickness=None, first_guess_t_bias=-2, t_bias_max_step_length=2, maxiter=30, store_model_geometry=True, store_fl_diagnostics=None, @@ -966,26 +979,27 @@ def dynamic_melt_f_run_with_dynamic_spinup( with the provided parameter climate_input_filesuffix and during the period of spinup_start_yr until rgi_year (e.g. 1979 - 2000). minimise_for : str - The variable we want to match at yr_rgi. Options are 'area' or 'volume'. + The variable we want to match at yr_target. Options are 'area' or + 'volume'. Default is 'area'. climate_input_filesuffix : str filesuffix for the input climate file Default is '' spinup_period : int The period how long the spinup should run. Start date of historical run - is defined "yr_rgi - spinup_period". Minimum allowed value is defined - with 'min_spinup_period'. If the provided climate data starts at year - later than (yr_rgi - spinup_period) the spinup_period is set to - (yr_rgi - yr_climate_start). Caution if spinup_start_yr is set the - spinup_period is ignored. + is defined "yr_target - spinup_period". Minimum allowed value is + defined with 'min_spinup_period'. If the provided climate data starts + at year later than (yr_target - spinup_period) the spinup_period is set + to (yr_target - yr_climate_start). Caution if spinup_start_yr is set + the spinup_period is ignored. Default is 20 min_spinup_period : int If the dynamic spinup function fails with the initial 'spinup_period' a shorter period is tried. Here you can define the minimum period to try. Default is 10 - yr_rgi : int or None - The rgi date, at which we want to match area or volume. + yr_target : int or None + The target year at which we want to match area or volume. If None, gdir.rgi_date + 1 is used (the default). Default is None precision_percent : float @@ -1049,8 +1063,8 @@ def dynamic_melt_f_run_with_dynamic_spinup( spinup_start_yr_max : int or None Possibility to provide a maximum year where the dynamic spinup must start from at least. If set, this overrides the min_spinup_period if - yr_rgi - spinup_start_yr_max > min_spinup_period. If None it is set to - yr0_ref_mb. + yr_target - spinup_start_yr_max > min_spinup_period. If None it is set + to yr0_ref_mb. Default is None add_fixed_geometry_spinup : bool If True and the original spinup_period of the dynamical spinup must be @@ -1087,15 +1101,15 @@ def dynamic_melt_f_run_with_dynamic_spinup( # we are done with preparing the local_variables for the upcoming iterations return None - if yr_rgi is None: - yr_rgi = gdir.rgi_date + 1 # + 1 converted to hydro years - if min_spinup_period > yr_rgi - ys: - log.info('The RGI year is closer to ys as the minimum spinup ' + if yr_target is None: + yr_target = gdir.rgi_date + 1 # + 1 converted to hydro years + if min_spinup_period > yr_target - ys: + log.info('The target year is closer to ys as the minimum spinup ' 'period -> therefore the minimum spinup period is ' 'adapted and it is the only period which is tried by the ' 'dynamic spinup function!') - min_spinup_period = yr_rgi - ys - spinup_period = yr_rgi - ys + min_spinup_period = yr_target - ys + spinup_period = yr_target - ys if spinup_start_yr_max is None: spinup_start_yr_max = yr0_ref_mb @@ -1162,7 +1176,7 @@ def dynamic_melt_f_run_with_dynamic_spinup( spinup_period=spinup_period, spinup_start_yr=ys, spinup_start_yr_max=spinup_start_yr_max, - min_spinup_period=min_spinup_period, yr_rgi=yr_rgi, + min_spinup_period=min_spinup_period, yr_target=yr_target, precision_percent=precision_percent, precision_absolute=precision_absolute, min_ice_thickness=min_ice_thickness, @@ -1204,7 +1218,7 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback( evolution_model=None, minimise_for='area', mb_model_historical=None, mb_model_spinup=None, climate_input_filesuffix='', spinup_period=20, min_spinup_period=10, - yr_rgi=None, precision_percent=1, + yr_target=None, precision_percent=1, precision_absolute=1, min_ice_thickness=10, first_guess_t_bias=-2, t_bias_max_step_length=2, maxiter=30, store_model_geometry=True, store_fl_diagnostics=None, @@ -1250,25 +1264,26 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback( with the provided parameter climate_input_filesuffix and during the period of spinup_start_yr until rgi_year (e.g. 1979 - 2000). minimise_for : str - The variable we want to match at yr_rgi. Options are 'area' or 'volume'. + The variable we want to match at yr_target. Options are 'area' or + 'volume'. Default is 'area'. climate_input_filesuffix : str filesuffix for the input climate file Default is '' spinup_period : int The period how long the spinup should run. Start date of historical run - is defined "yr_rgi - spinup_period". Minimum allowed value is defined - with 'min_spinup_period'. If the provided climate data starts at year - later than (yr_rgi - spinup_period) the spinup_period is set to - (yr_rgi - yr_climate_start). Caution if spinup_start_yr is set the - spinup_period is ignored. + is defined "yr_target - spinup_period". Minimum allowed value is + defined with 'min_spinup_period'. If the provided climate data starts + at year later than (yr_target - spinup_period) the spinup_period is set + to (yr_target - yr_climate_start). Caution if spinup_start_yr is set + the spinup_period is ignored. Default is 20 min_spinup_period : int If the dynamic spinup function fails with the initial 'spinup_period' a shorter period is tried. Here you can define the minimum period to try. Default is 10 - yr_rgi : int or None + yr_target : int or None The rgi date, at which we want to match area or volume. If None, gdir.rgi_date + 1 is used (the default). Default is None @@ -1321,7 +1336,7 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback( spinup_start_yr_max : int or None Possibility to provide a maximum year where the dynamic spinup must start from at least. If set, this overrides the min_spinup_period if - yr_rgi - spinup_start_yr_max > min_spinup_period. + yr_target - spinup_start_yr_max > min_spinup_period. Default is None add_fixed_geometry_spinup : bool If True and the original spinup_period of the dynamical spinup must be @@ -1363,15 +1378,15 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback( os.remove(os.path.join(gdir.dir, 'model_flowlines_dyn_melt_f_calib.pkl')) - if yr_rgi is None: - yr_rgi = gdir.rgi_date + 1 # + 1 converted to hydro years - if min_spinup_period > yr_rgi - ys: + if yr_target is None: + yr_target = gdir.rgi_date + 1 # + 1 converted to hydro years + if min_spinup_period > yr_target - ys: log.info('The RGI year is closer to ys as the minimum spinup ' 'period -> therefore the minimum spinup period is ' 'adapted and it is the only period which is tried by the ' 'dynamic spinup function!') - min_spinup_period = yr_rgi - ys - spinup_period = yr_rgi - ys + min_spinup_period = yr_target - ys + spinup_period = yr_target - ys yr_clim_min = gdir.get_climate_info()['baseline_yr_0'] try: @@ -1388,7 +1403,7 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback( spinup_start_yr=ys, min_spinup_period=min_spinup_period, spinup_start_yr_max=spinup_start_yr_max, - yr_rgi=yr_rgi, + yr_target=yr_target, minimise_for=minimise_for, precision_percent=precision_percent, precision_absolute=precision_absolute, @@ -1449,7 +1464,7 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback( def dynamic_melt_f_run( gdir, melt_f, yr0_ref_mb, yr1_ref_mb, fls_init, ys, ye, output_filesuffix='', evolution_model=None, - local_variables=None, set_local_variables=False, yr_rgi=None, + local_variables=None, set_local_variables=False, yr_target=None, **kwargs): """ This function is one option for a 'run_function' for the @@ -1487,8 +1502,8 @@ def dynamic_melt_f_run( set_local_variables : bool Not needed in this function. Only here to be confirm with the use of this function in 'run_dynamic_melt_f_calibration'. - yr_rgi : int or None - The rgi year of the gdir. + yr_target : int or None + The target year for a potential dynamic spinup (not needed here). Default is None kwargs : dict kwargs to pass to the evolution_model instance @@ -1538,7 +1553,7 @@ def dynamic_melt_f_run( def dynamic_melt_f_run_fallback( gdir, melt_f, fls_init, ys, ye, local_variables, output_filesuffix='', - evolution_model=None, yr_rgi=None, **kwargs): + evolution_model=None, yr_target=None, **kwargs): """ This is the fallback function corresponding to the function 'dynamic_melt_f_run', which are provided to @@ -1567,8 +1582,8 @@ def dynamic_melt_f_run_fallback( evolution_model : :class:oggm.core.FlowlineModel which evolution model to use. Default: cfg.PARAMS['evolution_model'] Not all models work in all circumstances! - yr_rgi : int or None - The rgi year of the gdir. + yr_target : int or None + The target year for a potential dynamic spinup (not needed here). Default is None kwargs : dict kwargs to pass to the evolution_model instance @@ -1608,7 +1623,7 @@ def run_dynamic_melt_f_calibration( ref_period='', melt_f_min=None, melt_f_max=None, melt_f_max_step_length_minimum=0.1, maxiter=20, ignore_errors=False, output_filesuffix='_dynamic_melt_f', - ys=None, ye=None, + ys=None, ye=None, yr_target=None, run_function=dynamic_melt_f_run_with_dynamic_spinup, kwargs_run_function=None, fallback_function=dynamic_melt_f_run_with_dynamic_spinup_fallback, @@ -1699,6 +1714,11 @@ def run_dynamic_melt_f_calibration( The end year of the conducted run. If None the last year of the provided climate file. Default is None + yr_target : int or None + The target year for a potential dynamic spinup (see run_dynamic_spinup + function for more info). + If None, gdir.rgi_date + 1 is used (the default). + Default is None run_function : function This function defines how a new defined melt_f is used to conduct the next model run. This function must contain the arguments 'gdir', @@ -1800,17 +1820,18 @@ def run_dynamic_melt_f_calibration( raise RuntimeError('The provided ys is larger than the start year of ' 'the given geodetic_mb_period!') - yr_rgi = gdir.rgi_date + 1 # + 1 converted to hydro years - if yr_rgi < ys: + if yr_target is None: + yr_target = gdir.rgi_date + 1 # + 1 converted to hydro years + if yr_target < ys: if ignore_errors: log.info('The rgi year is smaller than the provided start year ' 'ys -> setting the rgi year to ys to continue!') - yr_rgi = ys + yr_target = ys else: raise RuntimeError('The rgi year is smaller than the provided ' 'start year ys!') - kwargs_run_function['yr_rgi'] = yr_rgi - kwargs_fallback_function['yr_rgi'] = yr_rgi + kwargs_run_function['yr_target'] = yr_target + kwargs_fallback_function['yr_target'] = yr_target # get initial flowlines from which we want to start from if init_model_filesuffix is not None: diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index 48ecaaf8c..d6e23b602 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -3169,7 +3169,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): for store_model_evolution in [True, False]: model_dynamic_spinup = run_dynamic_spinup( hef_gdir, - yr_rgi=yr_rgi, + yr_target=yr_rgi, minimise_for=minimise_for, precision_percent=precision_percent, precision_absolute=precision_absolute, @@ -3229,13 +3229,48 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): assert fmod.last_yr == yr_rgi assert len(model_dynamic_spinup.fls) == len(fmod.fls) + # test user provided target year and value + yr_target = 2000 + if minimise_for == 'area': + ref_value = 8.5 + elif minimise_for == 'volume': + ref_value = 0.6 + model_dynamic_spinup_target_yr = run_dynamic_spinup( + hef_gdir, + yr_target=yr_target, + target_value=ref_value, + minimise_for=minimise_for, + precision_percent=precision_percent, + precision_absolute=precision_absolute, + min_ice_thickness=min_ice_thickness, + output_filesuffix='_dynamic_spinup', + store_model_evolution=store_model_evolution) + + # check if resulting model match wanted value with prescribed precision + if var_name == 'area_km2': + model_value = np.sum( + [np.sum(fl.bin_area_m2[fl.thick > min_ice_thickness]) + for fl in model_dynamic_spinup_target_yr.fls]) * 1e-6 + elif var_name == 'volume_km3': + model_value = np.sum( + [np.sum((fl.section * fl.dx_meter)[fl.thick > min_ice_thickness]) + for fl in model_dynamic_spinup_target_yr.fls]) * 1e-9 + else: + raise NotImplementedError(f'{var_name}') + assert np.isclose(model_value, ref_value, + rtol=precision_percent / 100, atol=0) + assert np.isclose(model_value, ref_value, + rtol=0, atol=precision_absolute) + assert model_dynamic_spinup_target_yr.yr == yr_target + assert len(model_dynamic_spinup_target_yr.fls) == len(fls) + # test if spinup_start_yr is handled correctly and overrides the spinup_period spinup_start_yr = yr_rgi - 20 model_dynamic_spinup_ys = run_dynamic_spinup( hef_gdir, spinup_period=40, spinup_start_yr=spinup_start_yr, - yr_rgi=yr_rgi, + yr_target=yr_rgi, minimise_for=minimise_for, precision_percent=precision_percent, precision_absolute=precision_absolute, @@ -3280,7 +3315,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): run_dynamic_spinup( hef_gdir, minimise_for=minimise_for, - yr_rgi=2002, + yr_target=2002, ye=2002, ignore_errors=ignore_errors, spinup_period=10, @@ -3299,7 +3334,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): model = run_dynamic_spinup( hef_gdir, minimise_for=minimise_for, - yr_rgi=2002, + yr_target=2002, ye=2002, ignore_errors=ignore_errors, maxiter=2, @@ -3324,7 +3359,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): hef_gdir, spinup_period=40, spinup_start_yr=spinup_start_yr, - yr_rgi=yr_rgi, + yr_target=yr_rgi, ye=ye, return_t_bias_best=True, minimise_for=minimise_for, @@ -3357,7 +3392,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): spinup_period=5, spinup_start_yr=None, spinup_start_yr_max=1990, - yr_rgi=yr_rgi, + yr_target=yr_rgi, minimise_for=minimise_for, precision_percent=precision_percent, precision_absolute=precision_absolute, @@ -3384,13 +3419,13 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): spinup_start_yr_max=yr_rgi - 10, spinup_start_yr=yr_rgi - 5) - # test that provided ye is larger than yr_rgi + # test that provided ye is larger than yr_target with pytest.raises(RuntimeError, match='The provided end year *'): run_dynamic_spinup( hef_gdir, minimise_for=minimise_for, - yr_rgi=yr_rgi, + yr_target=yr_rgi, ye=yr_rgi - 1) # test if provided model geometry works and some other principle @@ -3406,7 +3441,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): minimise_for=minimise_for, init_model_filesuffix='_one_yr', init_model_yr=yr_rgi - 1, - yr_rgi=yr_rgi, + yr_target=yr_rgi, store_model_geometry=False) # test that error is raised if mb_elev_feedback not annual @@ -3436,7 +3471,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): precision_percent=0.00012, minimise_for=minimise_for, output_filesuffix='_without_fixed_spinup', - yr_rgi=yr_rgi, + yr_target=yr_rgi, add_fixed_geometry_spinup=False) run_without_fixed_spinup = utils.compile_run_output( hef_gdir, input_filesuffix='_without_fixed_spinup', path=False) @@ -3448,7 +3483,7 @@ def test_run_dynamic_spinup(self, hef_gdir, minimise_for): precision_percent=0.00012, minimise_for=minimise_for, output_filesuffix='_with_fixed_spinup', - yr_rgi=yr_rgi, + yr_target=yr_rgi, add_fixed_geometry_spinup=True) run_with_fixed_spinup = utils.compile_run_output( hef_gdir, input_filesuffix='_with_fixed_spinup', path=False) From 6d6e8f31146694d571575a5a17450c69f1ec26c6 Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 26 Jun 2023 10:57:13 +0200 Subject: [PATCH 2/3] added target year to diagnostics --- oggm/core/dynamic_spinup.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/oggm/core/dynamic_spinup.py b/oggm/core/dynamic_spinup.py index a2f2bb632..bbf05d308 100644 --- a/oggm/core/dynamic_spinup.py +++ b/oggm/core/dynamic_spinup.py @@ -864,6 +864,8 @@ def get_mismatch(t_bias): # also save some other stuff gdir.add_to_diagnostics('temp_bias_dynamic_spinup', float(final_t_bias_guess[-1])) + gdir.add_to_diagnostics('dynamic_spinup_target_year', + int(yr_target)) gdir.add_to_diagnostics('dynamic_spinup_period', int(spinup_period)) gdir.add_to_diagnostics('dynamic_spinup_forward_model_iterations', From 1a47b57ed297d17b32670a8eba30d1b398ae4f09 Mon Sep 17 00:00:00 2001 From: Patrick Date: Mon, 26 Jun 2023 13:34:11 +0200 Subject: [PATCH 3/3] added whats-new --- docs/whats-new.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/whats-new.rst b/docs/whats-new.rst index f87fa3806..e1cffba0c 100644 --- a/docs/whats-new.rst +++ b/docs/whats-new.rst @@ -31,15 +31,15 @@ Enhancements functionality is currently in the sandbox and needs to be documented. By `Anouk Vlug `_ and `Fabien Maussion `_ - -Enhancements -~~~~~~~~~~~~ - - Added three new flowline diagnostic variables: thickness change in one year (``dhdt``), forcing climatic mass balance (``climatic_mb``) and flux divergence (``flux_divergence``). All variables are in units meter of ice per year (:pull:`1595`). By `Patrick Schmitt `_ +- Added more flexibility to ``run_dynamic_spinup``. Users can now specify a target + year and a desired value to match. The default is still the same, matching area + or volume at the RGI date (:pull:`1600`). + By `Patrick Schmitt `_ v1.6.0 (March 10, 2023)