From 1a97f3b2996a42bacd8dc3617f1e6bd7a3d7af92 Mon Sep 17 00:00:00 2001 From: Patrick Schmitt <50843444+pat-schmitt@users.noreply.github.com> Date: Tue, 4 Oct 2022 14:18:51 +0200 Subject: [PATCH] Introduce scaling factor for err_dmdtda and suppress logging, both in dynamic mu star calibration (#1488) * added err_dmdtda_scaling_factor in dynamic mu star calibration to reduce individual uncertainty * suppress logging of intermediate steps in dynamic mu star calibration * switched off logging also in fallback function of dynamic mu star calibration --- oggm/cli/prepro_levels.py | 17 ++++++-- oggm/core/dynamic_spinup.py | 67 ++++++++++++++++++++++-------- oggm/tests/test_models.py | 83 +++++++++++++++++++++++++++++++++++++ oggm/utils/_workflow.py | 3 +- oggm/workflow.py | 46 +++++++++++++------- 5 files changed, 180 insertions(+), 36 deletions(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 22ecf4aff..6383b7264 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -84,7 +84,8 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None, add_consensus=False, start_level=None, start_base_url=None, max_level=5, ref_tstars_base_url='', logging_level='WORKFLOW', disable_dl_verify=False, - dynamic_spinup=False, dynamic_spinup_start_year=1979, + dynamic_spinup=False, err_dmdtda_scaling_factor=1, + dynamic_spinup_start_year=1979, continue_on_error=True): """Generate the preprocessed OGGM glacier directories for this OGGM version @@ -154,8 +155,11 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None, a dict of parameters to override. disable_dl_verify : bool disable the hash verification of OGGM downloads - dynamic_spinup: str - include a dynamic spinup matching 'area' OR 'volume' at the RGI-date + dynamic_spinup : str + include a dynamic spinup matching 'area/dmdtda' OR 'volume/dmdtda' at + the RGI-date + err_dmdtda_scaling_factor : float + scaling factor to reduce individual geodetic mass balance uncertainty dynamic_spinup_start_year : int if dynamic_spinup is set, define the starting year for the simulation. The default is 1979, unless the climate data starts later. @@ -645,6 +649,7 @@ def _time_log(): workflow.execute_entity_task( tasks.run_dynamic_mu_star_calibration, gdirs, + err_dmdtda_scaling_factor=err_dmdtda_scaling_factor, ys=dynamic_spinup_start_year, ye=ye, max_mu_star=used_max_mu_star, kwargs_run_function={'evolution_model': evolution_model, @@ -835,6 +840,11 @@ def parse_args(args): "the RGI-date, AND mass-change from Hugonnet " "in the period 2000-2019 (dynamic mu* " "calibration).") + parser.add_argument('--err-dmdtda-scaling-factor', type=float, default=1, + help="scaling factor to account for correlated " + "uncertainties of geodetic mass balance " + "observations when looking at regional scale. " + "Should be smaller or equal 1.") parser.add_argument('--dynamic-spinup-start-year', type=int, default=1979, help="if --dynamic-spinup is set, define the starting" "year for the simulation. The default is 1979, " @@ -894,6 +904,7 @@ def parse_args(args): ref_tstars_base_url=args.ref_tstars_base_url, evolution_model=args.evolution_model, dynamic_spinup=dynamic_spinup, + err_dmdtda_scaling_factor=args.err_dmdtda_scaling_factor, dynamic_spinup_start_year=args.dynamic_spinup_start_year, ) diff --git a/oggm/core/dynamic_spinup.py b/oggm/core/dynamic_spinup.py index 693708f39..a9120dd26 100644 --- a/oggm/core/dynamic_spinup.py +++ b/oggm/core/dynamic_spinup.py @@ -1084,17 +1084,22 @@ def dynamic_mu_star_run_with_dynamic_spinup( define_new_mu_star_in_gdir(gdir, mu_star) if do_inversion: - apparent_mb_from_any_mb(gdir) - # do inversion with A calibration to current volume - calibrate_inversion_from_consensus( - [gdir], apply_fs_on_mismatch=True, error_on_mismatch=False, - filter_inversion_output=True, - volume_m3_reference=local_variables['vol_m3_ref']) + with utils.DisableLogger(): + apparent_mb_from_any_mb(gdir, + add_to_log_file=False, # dont write to log + ) + # do inversion with A calibration to current volume + calibrate_inversion_from_consensus( + [gdir], apply_fs_on_mismatch=True, error_on_mismatch=False, + filter_inversion_output=True, + volume_m3_reference=local_variables['vol_m3_ref'], + add_to_log_file=False) # this is used to keep the original model_flowline unchanged (-> to be able # to conduct different dynamic calibration runs in the same gdir) model_flowline_filesuffix = '_dyn_mu_calib' - init_present_time_glacier(gdir, filesuffix=model_flowline_filesuffix) + init_present_time_glacier(gdir, filesuffix=model_flowline_filesuffix, + add_to_log_file=False) # Now do a dynamic spinup to match area # do not ignore errors in dynamic spinup, so all 'bad' files are @@ -1103,6 +1108,7 @@ def dynamic_mu_star_run_with_dynamic_spinup( model, last_best_t_bias = run_dynamic_spinup( gdir, continue_on_error=False, # force to raise an error in @entity_task + add_to_log_file=False, # dont write to log file in @entity_task init_model_fls=fls_init, climate_input_filesuffix=climate_input_filesuffix, evolution_model=evolution_model, @@ -1136,8 +1142,9 @@ def dynamic_mu_star_run_with_dynamic_spinup( raise RuntimeError(f'Dynamic spinup raised error! (Message: {e})') # calculate dmdtda from previous simulation here - ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix, - path=False) + with utils.DisableLogger(): + ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix, + path=False) dmdtda_mdl = ((ds.volume.loc[yr1_ref_mb].values - ds.volume.loc[yr0_ref_mb].values) / gdir.rgi_area_m2 / @@ -1283,11 +1290,14 @@ def dynamic_mu_star_run_with_dynamic_spinup_fallback( if mu_star != gdir.read_json('local_mustar')['mu_star_glacierwide']: define_new_mu_star_in_gdir(gdir, mu_star) if do_inversion: - apparent_mb_from_any_mb(gdir) - calibrate_inversion_from_consensus( - [gdir], apply_fs_on_mismatch=True, error_on_mismatch=False, - filter_inversion_output=True, - volume_m3_reference=local_variables['vol_m3_ref']) + with utils.DisableLogger(): + apparent_mb_from_any_mb(gdir, + add_to_log_file=False) + calibrate_inversion_from_consensus( + [gdir], apply_fs_on_mismatch=True, error_on_mismatch=False, + filter_inversion_output=True, + volume_m3_reference=local_variables['vol_m3_ref'], + add_to_log_file=False) if os.path.isfile(os.path.join(gdir.dir, 'model_flowlines_dyn_mu_calib.pkl')): os.remove(os.path.join(gdir.dir, @@ -1309,6 +1319,7 @@ def dynamic_mu_star_run_with_dynamic_spinup_fallback( model_end = run_dynamic_spinup( gdir, continue_on_error=False, # force to raise an error in @entity_task + add_to_log_file=False, init_model_fls=fls_init, climate_input_filesuffix=climate_input_filesuffix, evolution_model=evolution_model, @@ -1341,7 +1352,9 @@ def dynamic_mu_star_run_with_dynamic_spinup_fallback( 'try is to conduct a run until ye without a dynamic ' 'spinup.') model_end = run_from_climate_data( - gdir, min_ys=yr_clim_min, ye=ye, + gdir, + add_to_log_file=False, + min_ys=yr_clim_min, ye=ye, output_filesuffix=output_filesuffix, climate_input_filesuffix=climate_input_filesuffix, store_model_geometry=store_model_geometry, @@ -1438,6 +1451,7 @@ def dynamic_mu_star_run( model = run_from_climate_data(gdir, # force to raise an error in @entity_task continue_on_error=False, + add_to_log_file=False, ys=ys, ye=ye, output_filesuffix=output_filesuffix, init_model_fls=fls_init, @@ -1448,8 +1462,9 @@ def dynamic_mu_star_run( f'(Message: {e})') # calculate dmdtda from previous simulation here - ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix, - path=False) + with utils.DisableLogger(): + ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix, + path=False) dmdtda_mdl = ((ds.volume.loc[yr1_ref_mb].values - ds.volume.loc[yr0_ref_mb].values) / gdir.rgi_area_m2 / @@ -1509,6 +1524,7 @@ def dynamic_mu_star_run_fallback( model = run_from_climate_data(gdir, # force to raise an error in @entity_task continue_on_error=False, + add_to_log_file=False, ys=ys, ye=ye, output_filesuffix=output_filesuffix, init_model_fls=fls_init, @@ -1524,7 +1540,7 @@ def dynamic_mu_star_run_fallback( @entity_task(log, writes=['inversion_flowlines']) def run_dynamic_mu_star_calibration( - gdir, ref_dmdtda=None, err_ref_dmdtda=None, + gdir, ref_dmdtda=None, err_ref_dmdtda=None, err_dmdtda_scaling_factor=1, ref_period='', ignore_hydro_months=False, min_mu_star=None, max_mu_star=None, mu_star_max_step_length=5, maxiter=20, ignore_errors=False, output_filesuffix='_dynamic_mu_star', @@ -1563,6 +1579,20 @@ def run_dynamic_mu_star_calibration( yr-1). Must always be a positive number. If None the data from Hugonett 2021 is used. Default is None + err_dmdtda_scaling_factor : float + The error of the geodetic mass balance is multiplied by this factor. + When looking at more glaciers you should set this factor smaller than + 1 (Default), but the smaller this factor the more glaciers will fail + during calibration. The factor is only used if ref_dmdtda = None and + err_ref_dmdtda = None. + The idea is that we reduce the uncertainty of individual observations + to count for correlated uncertainties when looking at regional or + global scales. If err_scaling_factor is 1 (Default) and you look at the + results of more than one glacier this equals that all errors are + uncorrelated. Therefore the result will be outside the uncertainty + boundaries given in Hugonett 2021 e.g. for the global estimate, because + some correlation of the individual errors is assumed during aggregation + of glaciers to regions (for more details see paper Hugonett 2021). ref_period : str If ref_dmdtda is None one of '2000-01-01_2010-01-01', '2010-01-01_2020-01-01', '2000-01-01_2020-01-01'. If ref_dmdtda is @@ -1697,6 +1727,7 @@ def run_dynamic_mu_star_calibration( err_ref_dmdtda = float(df_ref_dmdtda.loc[df_ref_dmdtda['period'] == ref_period]['err_dmdtda']) err_ref_dmdtda *= 1000 # kg m-2 yr-1 + err_ref_dmdtda *= err_dmdtda_scaling_factor if err_ref_dmdtda <= 0: raise RuntimeError('The provided error for the geodetic mass-balance ' diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index 96e5fa20f..239213a8b 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -3835,6 +3835,57 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self, for fl_prev, fl_now in zip(fls, gdir.read_pickle('model_flowlines'))]) + # test err_dmdtda_scaling_factor (not working for volume with inversion) + if not (do_inversion and minimise_for == 'volume'): + err_dmdtda_scaling_factor = 0.2 + run_dynamic_mu_star_calibration( + gdir, max_mu_star=1000., + err_dmdtda_scaling_factor=err_dmdtda_scaling_factor, + run_function=dynamic_mu_star_run_with_dynamic_spinup, + kwargs_run_function={'minimise_for': minimise_for, + 'precision_percent': precision_percent, + 'precision_absolute': precision_absolute, + 'do_inversion': do_inversion}, + fallback_function=dynamic_mu_star_run_with_dynamic_spinup_fallback, + kwargs_fallback_function={'minimise_for': minimise_for, + 'precision_percent': precision_percent, + 'precision_absolute': precision_absolute, + 'do_inversion': do_inversion}, + output_filesuffix='_dyn_mu_calib_err_scaling', + ys=1979, ye=ye) + + # check that we are matching all desired ref values + ds = utils.compile_run_output( + gdir, input_filesuffix='_dyn_mu_calib_err_scaling', + path=False) + if minimise_for == 'volume': + assert np.isclose(ds.loc[{'time': yr_rgi}].volume.values * 1e-9, + ref_value_dynamic_spinup, + rtol=precision_percent / 100, + atol=precision_absolute) + elif minimise_for == 'area': + assert np.isclose(ds.loc[{'time': yr_rgi}].area.values * 1e-6, + ref_value_dynamic_spinup, + rtol=precision_percent / 100, + atol=precision_absolute) + dmdtda_mdl_scale = ((ds.volume.loc[yr1_ref_dmdtda].values - + ds.volume.loc[yr0_ref_dmdtda].values) / + gdir.rgi_area_m2 / + (yr1_ref_dmdtda - yr0_ref_dmdtda) * + cfg.PARAMS['ice_density']) + assert np.isclose(dmdtda_mdl_scale, ref_dmdtda, + rtol=np.abs(err_ref_dmdtda * + err_dmdtda_scaling_factor / ref_dmdtda)) + # check that calibration without scaling factor is outside adapted + # uncertainty (use previous result without err scaling factor), this + # tests if the scaling factor has any effect + assert not np.isclose(dmdtda_mdl, ref_dmdtda, + rtol=np.abs(err_ref_dmdtda * + err_dmdtda_scaling_factor / + ref_dmdtda)) + assert gdir.get_diagnostics()['used_spinup_option'] == \ + 'dynamic mu_star calibration (full success)' + # test that error is raised if user provides flowlines but want to # include inversion during dynamic mu calibration if do_inversion: @@ -4120,6 +4171,38 @@ def test_run_dynamic_mu_star_calibration_without_dynamic_spinup(self): assert gdir.get_diagnostics()['used_spinup_option'] == \ 'dynamic mu_star calibration (full success)' + # test err_dmdtda_scaling_factor + err_dmdtda_scaling_factor = 0.01 + run_dynamic_mu_star_calibration( + gdir, max_mu_star=1000., + err_dmdtda_scaling_factor=err_dmdtda_scaling_factor, + run_function=dynamic_mu_star_run, + fallback_function=dynamic_mu_star_run_fallback, + output_filesuffix='_dyn_mu_calib_err_scaling', + ys=1979, ye=ye) + # check that we are matching all desired ref values + ds = utils.compile_run_output( + gdir, input_filesuffix='_dyn_mu_calib_err_scaling', + path=False) + + dmdtda_mdl_scale = ((ds.volume.loc[yr1_ref_dmdtda].values - + ds.volume.loc[yr0_ref_dmdtda].values) / + gdir.rgi_area_m2 / + (yr1_ref_dmdtda - yr0_ref_dmdtda) * + cfg.PARAMS['ice_density']) + assert np.isclose(dmdtda_mdl_scale, ref_dmdtda, + rtol=np.abs(err_ref_dmdtda * + err_dmdtda_scaling_factor / ref_dmdtda)) + # check that calibration without scaling factor is outside adapted + # uncertainty (use previous result without err scaling factor), this + # tests if the scaling factor has any effect + assert not np.isclose(dmdtda_mdl, ref_dmdtda, + rtol=np.abs(err_ref_dmdtda * + err_dmdtda_scaling_factor / + ref_dmdtda)) + assert gdir.get_diagnostics()['used_spinup_option'] == \ + 'dynamic mu_star calibration (full success)' + # tests for user provided dmdtda (always reset gdir before each test) gdir = workflow.init_glacier_directories( ['RGI60-11.00897'], # Hintereisferner diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py index 01a1c5c5b..c3defc123 100644 --- a/oggm/utils/_workflow.py +++ b/oggm/utils/_workflow.py @@ -493,7 +493,8 @@ def _entity_task(gdir, *, reset=None, print_log=True, if cfg.PARAMS['task_timeout'] > 0: signal.alarm(0) if task_name != 'gdir_to_tar': - gdir.log(task_name, task_time=ex_t) + if add_to_log_file: + gdir.log(task_name, task_time=ex_t) except Exception as err: # Something happened out = None diff --git a/oggm/workflow.py b/oggm/workflow.py index f98ea9fc8..339e25fdd 100644 --- a/oggm/workflow.py +++ b/oggm/workflow.py @@ -502,7 +502,8 @@ def climate_tasks(gdirs, base_url=None): @global_task(log) -def inversion_tasks(gdirs, glen_a=None, fs=None, filter_inversion_output=True): +def inversion_tasks(gdirs, glen_a=None, fs=None, filter_inversion_output=True, + add_to_log_file=True): """Run all ice thickness inversion tasks on a list of glaciers. Quite useful to deal with calving glaciers as well. @@ -511,6 +512,8 @@ def inversion_tasks(gdirs, glen_a=None, fs=None, filter_inversion_output=True): ---------- gdirs : list of :py:class:`oggm.GlacierDirectory` objects the glacier directories to process + add_to_log_file : bool + if the called entity tasks should write into log of gdir. Default True """ if cfg.PARAMS['use_kcalving_for_inversion']: @@ -528,21 +531,28 @@ def inversion_tasks(gdirs, glen_a=None, fs=None, filter_inversion_output=True): len(gdirs_nc))) if gdirs_nc: - execute_entity_task(tasks.prepare_for_inversion, gdirs_nc) + execute_entity_task(tasks.prepare_for_inversion, gdirs_nc, + add_to_log_file=add_to_log_file) execute_entity_task(tasks.mass_conservation_inversion, gdirs_nc, - glen_a=glen_a, fs=fs) + glen_a=glen_a, fs=fs, + add_to_log_file=add_to_log_file) if filter_inversion_output: - execute_entity_task(tasks.filter_inversion_output, gdirs_nc) + execute_entity_task(tasks.filter_inversion_output, gdirs_nc, + add_to_log_file=add_to_log_file) if gdirs_c: execute_entity_task(tasks.find_inversion_calving, gdirs_c, - glen_a=glen_a, fs=fs) + glen_a=glen_a, fs=fs, + add_to_log_file=add_to_log_file) else: - execute_entity_task(tasks.prepare_for_inversion, gdirs) + execute_entity_task(tasks.prepare_for_inversion, gdirs, + add_to_log_file=add_to_log_file) execute_entity_task(tasks.mass_conservation_inversion, gdirs, - glen_a=glen_a, fs=fs) + glen_a=glen_a, fs=fs, + add_to_log_file=add_to_log_file) if filter_inversion_output: - execute_entity_task(tasks.filter_inversion_output, gdirs) + execute_entity_task(tasks.filter_inversion_output, gdirs, + add_to_log_file=add_to_log_file) @global_task(log) @@ -551,7 +561,8 @@ def calibrate_inversion_from_consensus(gdirs, ignore_missing=True, apply_fs_on_mismatch=False, error_on_mismatch=True, filter_inversion_output=True, - volume_m3_reference=None): + volume_m3_reference=None, + add_to_log_file=True): """Fit the total volume of the glaciers to the 2019 consensus estimate. This method finds the "best Glen A" to match all glaciers in gdirs with @@ -580,6 +591,8 @@ def calibrate_inversion_from_consensus(gdirs, ignore_missing=True, output (needs the downstream lines to work). volume_m3_reference : float Option to give an own total glacier volume to match to + add_to_log_file : bool + if the called entity tasks should write into log of gdir. Default True Returns ------- @@ -606,9 +619,11 @@ def calibrate_inversion_from_consensus(gdirs, ignore_missing=True, def compute_vol(x): inversion_tasks(gdirs, glen_a=x*def_a, fs=fs, - filter_inversion_output=filter_inversion_output) + filter_inversion_output=filter_inversion_output, + add_to_log_file=add_to_log_file) odf = df.copy() - odf['oggm'] = execute_entity_task(tasks.get_inversion_volume, gdirs) + odf['oggm'] = execute_entity_task(tasks.get_inversion_volume, gdirs, + add_to_log_file=add_to_log_file) # if the user provides a glacier volume all glaciers are considered, # dropna() below exclude glaciers where no ITMIX volume is available if volume_m3_reference is None: @@ -657,7 +672,8 @@ def to_minimize(x): fs=5.7e-20, a_bounds=a_bounds, apply_fs_on_mismatch=False, error_on_mismatch=error_on_mismatch, - volume_m3_reference=volume_m3_reference) + volume_m3_reference=volume_m3_reference, + add_to_log_file=add_to_log_file) if error_on_mismatch: raise ValueError(msg) @@ -669,8 +685,10 @@ def to_minimize(x): # Compute the final volume with the correct A inversion_tasks(gdirs, glen_a=out_fac*def_a, fs=fs, - filter_inversion_output=filter_inversion_output) - df['vol_oggm_m3'] = execute_entity_task(tasks.get_inversion_volume, gdirs) + filter_inversion_output=filter_inversion_output, + add_to_log_file=add_to_log_file) + df['vol_oggm_m3'] = execute_entity_task(tasks.get_inversion_volume, gdirs, + add_to_log_file=add_to_log_file) return df