From 12a025fb84e10e6707fb50d22f7f9626f5b1fadd Mon Sep 17 00:00:00 2001 From: Patrick Schmitt <50843444+pat-schmitt@users.noreply.github.com> Date: Fri, 1 Sep 2023 13:30:21 +0200 Subject: [PATCH] New sandbox experiments (#38) * add three different glacier states per sandbox geometry * typo * fix tests * adapted dashboard for different glacier states * added initial_flux to data_logger * added tests for glacier_state mb_models * small experiment adaptions --- agile1d/core/cost_function.py | 53 +- agile1d/core/data_logging.py | 5 +- agile1d/core/inversion.py | 14 +- agile1d/core/utils.py | 4 +- agile1d/sandbox/calculate_statistics.py | 58 +- .../create_glaciers_with_measurements.py | 378 ++++++++--- .../sandbox/define_idealized_experiment.py | 93 ++- agile1d/sandbox/graphics.py | 170 +++-- .../individual_experiment_dashboard.py | 35 +- agile1d/sandbox/run_idealized_experiment.py | 5 + agile1d/tests/test_sandbox.py | 618 ++++++++++-------- 11 files changed, 930 insertions(+), 503 deletions(-) diff --git a/agile1d/core/cost_function.py b/agile1d/core/cost_function.py index e8c2c6e..3caf57d 100644 --- a/agile1d/core/cost_function.py +++ b/agile1d/core/cost_function.py @@ -138,7 +138,6 @@ def get_indices_for_unknown_parameters(data_logger): def define_scaling_terms(data_logger): - # define scaling terms for observation stuff if 'scale' in data_logger.obs_scaling_parameters.keys(): observations = data_logger.observations @@ -260,10 +259,7 @@ def cost_fct(unknown_parameters, data_logger): mb_models = initialise_mb_models(unknown_parameters_descaled, data_logger) - if 'smoothed_flux' in data_logger.regularisation_terms.keys(): - initial_flux = dynamic_model(flowline).flux_stag[0] - else: - initial_flux = None + initial_flux = dynamic_model(flowline).flux_stag[0] if data_logger.spinup_type == 'height_shift_spinup': try: @@ -318,7 +314,7 @@ def cost_fct(unknown_parameters, data_logger): initial_flux, # for regularisation term 'smoothed_flux' unknown_parameters, # for distance from fg regularisation data_logger # for reg_parameters and observations - ) + ) # sum up cost function terms using cost lambda cost_lambda = torch.tensor(data_logger.cost_lambda, @@ -344,6 +340,9 @@ def cost_fct(unknown_parameters, data_logger): data_logger.save_data_in_datalogger('sfc_h_start', sfc_h_start) data_logger.save_data_in_datalogger('section_start', section_start) data_logger.save_data_in_datalogger('flowlines', detach_flowline(final_fl)) + # exclude last grid point because initial flux is defined on staggered grid + data_logger.save_data_in_datalogger('initial_flux', + initial_flux.detach().clone()[:-1]) data_logger.save_data_in_datalogger('costs', cost) data_logger.save_data_in_datalogger('grads', grad) data_logger.save_data_in_datalogger('c_terms', c_terms) @@ -625,20 +624,36 @@ def initialise_mb_models(unknown_parameters, # -1 because period defined as [y0 - halfsize, y0 + halfsize + 1] y0 = (y_start + y_end - 1) / 2 halfsize = (y_end - y_start - 1) / 2 - mb_models[mb_mdl_set] = {'mb_model': - ConstantMassBalanceTorch( - gdir, y0=y0, halfsize=halfsize, - torch_type=torch_type, device=device), - 'years': - mb_models_settings[mb_mdl_set]['years']} + mb_models[mb_mdl_set] = { + 'mb_model': ConstantMassBalanceTorch( + gdir, y0=y0, halfsize=halfsize, + torch_type=torch_type, device=device), + 'years': + mb_models_settings[mb_mdl_set]['years']} elif mb_models_settings[mb_mdl_set]['type'] == 'TIModel': - mb_models[mb_mdl_set] = {'mb_model': - MBModelTorchWrapper( - gdir=gdir, - mb_model=MonthlyTIModel(gdir)), - 'years': - mb_models_settings[mb_mdl_set]['years'] - } + mb_model_args = mb_models_settings[mb_mdl_set]['model_args'] + mb_models[mb_mdl_set] = { + 'mb_model': MBModelTorchWrapper( + gdir=gdir, + mb_model=MonthlyTIModel( + gdir, + **mb_model_args) + ), + 'years': + mb_models_settings[mb_mdl_set]['years'] + } + elif mb_models_settings[mb_mdl_set]['type'] == 'ConstantModel': + mb_model_args = mb_models_settings[mb_mdl_set]['model_args'] + mb_models[mb_mdl_set] = { + 'mb_model': MBModelTorchWrapper( + gdir=gdir, + mb_model=ConstantMassBalance( + gdir, + **mb_model_args) + ), + 'years': + mb_models_settings[mb_mdl_set]['years'] + } else: raise NotImplementedError("The MassBalance type " f"{mb_models_settings[mb_mdl_set]['type']} " diff --git a/agile1d/core/data_logging.py b/agile1d/core/data_logging.py index 58358e5..57845ea 100644 --- a/agile1d/core/data_logging.py +++ b/agile1d/core/data_logging.py @@ -140,6 +140,7 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename, self.time_needed = None self.grads = None self.flowlines = None + self.initial_flux = None self.sfc_h_start = None self.section_start = None self.observations_mdl = None @@ -157,7 +158,7 @@ def __init__(self, gdir, fls_init, inversion_input, climate_filename, self.first_guess = None self.filename = gdir.name + '_' + \ - inversion_input['experiment_description'] + inversion_input['experiment_description'] # create info Text for callback_fct TODO: think about showing the evolution of the c_terms self.info_text = ''' @@ -248,6 +249,7 @@ def filter_data_from_optimization(self): self.reg_terms = self.reg_terms[index] self.time_needed = self.squeeze_generic(self.time_needed[index]) self.flowlines = self.squeeze_generic(self.flowlines[index]) + self.initial_flux = self.squeeze_generic(self.initial_flux[index]) self.sfc_h_start = self.squeeze_generic(self.sfc_h_start[index]) self.section_start = self.squeeze_generic(self.section_start[index]) self.observations_mdl = self.squeeze_generic(self.observations_mdl[index]) @@ -269,6 +271,7 @@ def create_and_save_dataset(self): ds['costs'] = (['iteration'], self.costs) ds['grads'] = (['iteration', 'nr_unknown_parameters'], self.grads) ds['flowlines'] = (['iteration'], self.flowlines) + ds['initial_flux'] = (['iteration', 'x'], self.initial_flux) ds['sfc_h_start'] = (['iteration', 'x'], self.sfc_h_start) ds['section_start'] = (['iteration', 'x'], self.section_start) ds['observations_mdl'] = (['iteration'], self.observations_mdl) diff --git a/agile1d/core/inversion.py b/agile1d/core/inversion.py index 9c34765..04f21e1 100644 --- a/agile1d/core/inversion.py +++ b/agile1d/core/inversion.py @@ -60,9 +60,11 @@ def add_setting(): "MassBalanceModel must start at least " \ "one year before first given observation year!" \ "Default: {'MB': {'type': 'TIModel'," \ - " 'years': np.array([1980, 2020])}}" + " 'years': np.array([1980, 2020])," \ + " 'model_args': {}}}" _default = {'MB': {'type': 'TIModel', - 'years': np.array([1980, 2020])}} + 'years': np.array([1980, 2020]), + 'model_args': {}}} add_setting() @@ -103,10 +105,10 @@ def add_setting(): _key = "bed_h_bounds" _doc = "Define how large the boundaries for the bed_h are, in relation of " \ - "first guess thickness. (e.g. (0.2, 1.4) means the bed height can " \ + "first guess thickness. (e.g. (0.4, 1.6) means the bed height can " \ "be between 1.4*fg_thick and 0.2*fg_thick). " \ - "Default: (0.2, 1.4)" - _default = (0.2, 1.4) + "Default: (0.4, 1.6)" + _default = (0.4, 1.6) add_setting() _key = "max_deviation_surface_h" @@ -245,7 +247,7 @@ def add_setting(): "glacier at the inital state." \ "e.g. {'section':" \ " {'extra_grid_points': 10," \ - " 'limits': (0.75, 1.25)," \ + " 'limits': (0.6, 1.4)," \ " 'fg_years': 1" \ " }" \ " }" \ diff --git a/agile1d/core/utils.py b/agile1d/core/utils.py index 820dd42..0d928a8 100644 --- a/agile1d/core/utils.py +++ b/agile1d/core/utils.py @@ -144,7 +144,9 @@ def inversion_settings(self): def write_inversion_settings(self, control_vars=['bed_h'], - mb_models_settings={'MB1': {'type': 'constant', 'years': np.array([1950, 2016])}}, + mb_models_settings={'MB1': { + 'type': 'constant', + 'years': np.array([1950, 2016])}}, min_w0_m=10., observations=None, # {'Area', {'2010': np.array([23])}} reg_parameters=None, # [0, 0.1, 10] diff --git a/agile1d/sandbox/calculate_statistics.py b/agile1d/sandbox/calculate_statistics.py index 166a92b..9abf22d 100644 --- a/agile1d/sandbox/calculate_statistics.py +++ b/agile1d/sandbox/calculate_statistics.py @@ -30,7 +30,8 @@ def add_0d_stats(x, y): 'abs_diff': float(np.abs(x - y))} -def calculate_result_statistics(gdir, data_logger, print_statistic=False): +def calculate_result_statistics(gdir, glacier_state, data_logger, + print_statistic=False): """calculate some statistics of the result for analysis""" # open the dataset of the run to add our calculated statistics @@ -69,7 +70,8 @@ def calculate_result_statistics(gdir, data_logger, print_statistic=False): ds.unknown_parameters[-1][control_indices[control_var]].values controls_true = {} fls_true = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] for control_var in controls_mdl: if control_var in ['bed_h', 'area_bed_h']: controls_true[control_var] = \ @@ -81,7 +83,8 @@ def calculate_result_statistics(gdir, data_logger, print_statistic=False): controls_true[control_var] = controls_mdl[control_var] elif control_var in ['section']: fls_1980_true = gdir.read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] + filesuffix='_creation_spinup_' + f'{glacier_state}')[0] controls_true[control_var] = \ fls_1980_true.section[:len(controls_mdl['section'])] else: @@ -117,7 +120,8 @@ def calculate_result_statistics(gdir, data_logger, print_statistic=False): sfc_h_start = copy.deepcopy(data_logger.sfc_h_start[-1]) fls_start_mdl.surface_h = sfc_h_start fls_start_true = gdir.read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] + filesuffix='_creation_spinup_' + f'{glacier_state}')[0] past_state_stats = {} for var in ['thick', 'area_m2', 'volume_m3']: @@ -150,7 +154,8 @@ def get_volume(fl): with xr.open_dataset(fp) as ds_diag: past_evol_mdl = ds_diag.load() fp = gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_total_run') + filesuffix='_agile_true_total_run_' + f'{glacier_state}') with xr.open_dataset(fp) as ds_diag: past_evol_true = ds_diag.load() @@ -165,7 +170,8 @@ def get_volume(fl): # how well do we match today's glacier state ------------------------------ fls_end_mdl = copy.deepcopy(data_logger.flowlines[-1]) fls_end_true = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_end')[0] + filesuffix='_agile_true_end_' + f'{glacier_state}')[0] today_state_stats = {} for var in ['thick', 'area_m2', 'volume_m3']: @@ -196,7 +202,8 @@ def get_volume(fl): with xr.open_dataset(fp) as ds_diag: future_evol_mdl = ds_diag.load() fp = gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_future') + filesuffix='_agile_true_future_' + f'{glacier_state}') with xr.open_dataset(fp) as ds_diag: future_evol_true = ds_diag.load() @@ -218,7 +225,8 @@ def get_volume(fl): # Here print the statistics in comparision to the default run # open default statistics - fp_default = os.path.join(gdir.dir, 'default_oggm_statistics.pkl') + fp_default = os.path.join(gdir.dir, f'default_oggm_statistics_' + f'{glacier_state}.pkl') with open(fp_default, 'rb') as handle: default_stats = pickle.load(handle) @@ -271,7 +279,7 @@ def get_volume(fl): raise NotImplementedError(f'{stat_clean}') -def calculate_default_oggm_statistics(gdir): +def calculate_default_oggm_statistics(gdir, glacier_state): default_oggm_statistics = {} @@ -279,20 +287,24 @@ def calculate_default_oggm_statistics(gdir): # open the run files with xr.open_dataset( gdir.get_filepath('model_diagnostics', - filesuffix=f'_oggm_{reali}_past')) as ds: + filesuffix=f'_oggm_{reali}_past_' + f'{glacier_state}')) as ds: diag_past = ds.load() f = gdir.get_filepath('fl_diagnostics', - filesuffix=f'_oggm_{reali}_past') + filesuffix=f'_oggm_{reali}_past_' + f'{glacier_state}') with xr.open_dataset(f, group=f'fl_0') as ds: fl_diag_past = ds.load() with xr.open_dataset( gdir.get_filepath('model_diagnostics', - filesuffix=f'_oggm_{reali}_future')) as ds: + filesuffix=f'_oggm_{reali}_future_' + f'{glacier_state}')) as ds: diag_future = ds.load() # how well do we match the observations ------------------------------- obs_given = gdir.read_pickle('inversion_input', - filesuffix='_agile_measurements') + filesuffix='_agile_measurements_' + f'{glacier_state}') obs_stats = {} for obs_key in obs_given.keys(): obs_stats[obs_key] = {} @@ -374,9 +386,11 @@ def calculate_default_oggm_statistics(gdir): controls_mdl = {} controls_true = {} fls_mdl = gdir.read_pickle('model_flowlines', - filesuffix='_oggm_first_guess')[0] + filesuffix='_oggm_first_guess_' + f'{glacier_state}')[0] fls_true = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] for control_var in all_control_vars: if control_var in ['bed_h', 'area_bed_h']: @@ -405,7 +419,8 @@ def calculate_default_oggm_statistics(gdir): # how well do we match the past glacier state ----------------------------- fls_start_mdl = fl_diag_past.sel(time=fl_diag_past.time[0]) fls_start_true = gdir.read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] + filesuffix='_creation_spinup_' + f'{glacier_state}')[0] past_state_stats = {} for var in ['thick', 'area_m2', 'volume_m3']: @@ -435,7 +450,8 @@ def get_volume(fl): # how well do we match the past glacier evolution --------------------- past_evol_mdl = diag_past fp = gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_total_run') + filesuffix='_agile_true_total_run_' + f'{glacier_state}') with xr.open_dataset(fp) as ds_diag: past_evol_true = ds_diag.load() @@ -450,7 +466,8 @@ def get_volume(fl): # how well do we match today's glacier state -------------------------- fls_end_mdl = fl_diag_past.sel(time=fl_diag_past.time[-1]) fls_end_true = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_end')[0] + filesuffix='_agile_true_end_' + f'{glacier_state}')[0] today_state_stats = {} for var in ['thick', 'area_m2', 'volume_m3']: @@ -480,7 +497,8 @@ def get_volume(fl): # how well do we match the future glacier evolution ------------------- future_evol_mdl = diag_future fp = gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_future') + filesuffix='_agile_true_future_' + f'{glacier_state}') with xr.open_dataset(fp) as ds_diag: future_evol_true = ds_diag.load() @@ -494,7 +512,7 @@ def get_volume(fl): # save final default statistics as pickle out = os.path.join(gdir.dir, - 'default_oggm_statistics.pkl') + f'default_oggm_statistics_{glacier_state}.pkl') with open(out, 'wb') as handle: pickle.dump(default_oggm_statistics, handle, protocol=pickle.HIGHEST_PROTOCOL) diff --git a/agile1d/sandbox/create_glaciers_with_measurements.py b/agile1d/sandbox/create_glaciers_with_measurements.py index ac5f82e..5403b70 100644 --- a/agile1d/sandbox/create_glaciers_with_measurements.py +++ b/agile1d/sandbox/create_glaciers_with_measurements.py @@ -21,6 +21,7 @@ def create_idealized_experiments(all_glaciers, + glacier_states, prepro_border=None, from_prepro_level=None, base_url=None, @@ -51,6 +52,14 @@ def create_idealized_experiments(all_glaciers, if glacier not in experiment_glaciers.keys(): raise ValueError(f'{glacier} not defined for idealized experiments') + # check if glacier states are valid + if any(state not in ['equilibrium', 'retreating', 'advancing'] + for state in glacier_states): + raise NotImplementedError("Glacier states could only be one of " + "['equilibrium', 'retreating', 'advancing']" + ", but given was " + f"{glacier_states}.") + # create a dict which translate rgi_ids to glacier name rgi_id_to_name = dict([(experiment_glaciers[glacier]['rgi_id'], glacier) for glacier in all_glaciers]) @@ -116,29 +125,36 @@ def create_idealized_experiments(all_glaciers, for gdir in gdirs: gdir.rgi_date = 2000 - # create spinup glacier starting from consensus flowline - workflow.execute_entity_task(create_spinup_glacier, gdirs, - yr_run_start=yr_run_start) - - # create glacier with experiment measurements - workflow.execute_entity_task(evolve_glacier_and_create_measurements, gdirs, - used_mb_models=used_mb_models, - yr_dmdt_start=yr_dmdt_start, - yr_run_start=yr_run_start, - yr_run_end=yr_run_end, - gcm=gcm, - ssp=ssp) - - # now OGGM inversion from idealized glacier surface for first guess - workflow.execute_entity_task(oggm_inversion_for_first_guess, gdirs) - - # now add a second bed inversion - workflow.execute_entity_task(glabtop_inversion_for_first_guess, gdirs) - - # finally use oggm default initialisation for comparisons - workflow.execute_entity_task(oggm_default_initialisation, gdirs, - ys=yr_run_start, ye=yr_run_end, gcm=gcm, - ssp=ssp) + # create glacier with experiment measurements, for all states + for glacier_state in glacier_states: + # create spinup glacier starting from consensus flowline + workflow.execute_entity_task(create_spinup_glacier, gdirs, + glacier_state=glacier_state, + yr_run_start=yr_run_start) + + workflow.execute_entity_task(evolve_glacier_and_create_measurements, + gdirs, + glacier_state=glacier_state, + used_mb_models=used_mb_models, + yr_dmdt_start=yr_dmdt_start, + yr_run_start=yr_run_start, + yr_run_end=yr_run_end, + gcm=gcm, + ssp=ssp) + + # now OGGM inversion from idealized glacier surface for first guess + workflow.execute_entity_task(oggm_inversion_for_first_guess, + gdirs, glacier_state=glacier_state) + + # now add a second bed inversion + workflow.execute_entity_task(glabtop_inversion_for_first_guess, + gdirs, glacier_state=glacier_state) + + # finally use oggm default initialisation for comparisons + workflow.execute_entity_task(oggm_default_initialisation, + gdirs, glacier_state=glacier_state, + ys=yr_run_start, ye=yr_run_end, gcm=gcm, + ssp=ssp) print('Experiment creation finished') return gdirs @@ -234,8 +250,75 @@ def initialise_consensus_flowline(gdir): filesuffix='_consensus') +def get_equilibrium_year_and_halfsize(gdir): + return ({'RGI60-11.01450': 1993, # Aletsch + 'RGI60-16.02444': 1982, # Artesonraju + 'RGI60-14.06794': 1987, # Baltoro + 'RGI60-02.05098': 1957, # Peyto + }[gdir.rgi_id], + {'RGI60-11.01450': 5, # Aletsch + 'RGI60-16.02444': 1, # Artesonraju + 'RGI60-14.06794': 6, # Baltoro + 'RGI60-02.05098': 9, # Peyto + }[gdir.rgi_id]) + + +def get_spinup_mb_model(gdir, glacier_state, fls_spinup): + if glacier_state in ['retreating']: + t_bias_default = gdir.read_json('mb_calib')['temp_bias'] + if gdir.rgi_id == 'RGI60-16.02444': + t_bias_use = t_bias_default + 0 + elif gdir.rgi_id in ['RGI60-02.05098', 'RGI60-14.06794']: + t_bias_use = t_bias_default - 0 + elif gdir.rgi_id == 'RGI60-11.01450': + t_bias_use = t_bias_default + 1.1 + else: + raise NotImplementedError(f'{gdir.rgi_id}') + mb_model = MultipleFlowlineMassBalance(gdir, + fls=fls_spinup, + mb_model_class=RandomMassBalance, + seed=2, + filename='climate_historical', + y0=1950, + halfsize=30, + temp_bias=t_bias_use + ) + elif glacier_state in ['advancing']: + t_bias_default = gdir.read_json('mb_calib')['temp_bias'] + if gdir.rgi_id == 'RGI60-16.02444': + t_bias_use = t_bias_default + 0.4 + elif gdir.rgi_id in ['RGI60-02.05098']: + t_bias_use = t_bias_default - 0.5 + elif gdir.rgi_id in ['RGI60-14.06794']: + t_bias_use = t_bias_default - 0 + elif gdir.rgi_id == 'RGI60-11.01450': + t_bias_use = t_bias_default + 1.1 + else: + raise NotImplementedError(f'{gdir.rgi_id}') + mb_model = MultipleFlowlineMassBalance(gdir, + fls=fls_spinup, + mb_model_class=RandomMassBalance, + seed=2, + filename='climate_historical', + y0=1950, + halfsize=30, + temp_bias=t_bias_use) + elif glacier_state in ['equilibrium']: + y0, halfsize = get_equilibrium_year_and_halfsize(gdir) + mb_model = MultipleFlowlineMassBalance(gdir, + fls=fls_spinup, + mb_model_class=ConstantMassBalance, + filename='climate_historical', + y0=y0, + halfsize=halfsize) + else: + raise NotImplementedError(f'{glacier_state}') + + return mb_model + + @entity_task(log, writes=['model_flowlines']) -def create_spinup_glacier(gdir, yr_run_start): +def create_spinup_glacier(gdir, glacier_state, yr_run_start): """ Creates and saves the glacier state at yr_run_start (e.g. 1980) as 'model_flowlines_creation_spinup'. Method inspired from model creation of @@ -243,6 +326,8 @@ def create_spinup_glacier(gdir, yr_run_start): """ # define how many years the random climate spinup should be years_to_run_random = 60 + if glacier_state == 'equilibrium': + years_to_run_random = 120 # use flow-parameters from default oggm inversion diag = gdir.get_diagnostics() @@ -254,13 +339,7 @@ def create_spinup_glacier(gdir, yr_run_start): fls_spinup = gdir.read_pickle('model_flowlines', filesuffix='_consensus') - mb_random = MultipleFlowlineMassBalance(gdir, - fls=fls_spinup, - mb_model_class=RandomMassBalance, - seed=2, - filename='climate_historical', - y0=1950, - halfsize=30) + mb_random = get_spinup_mb_model(gdir, glacier_state, fls_spinup) model = SemiImplicitModel(fls_spinup, mb_random, @@ -269,24 +348,121 @@ def create_spinup_glacier(gdir, yr_run_start): model.run_until_and_store( yr_run_start, diag_path=gdir.get_filepath('model_diagnostics', - filesuffix='_agile_creation_spinup', + filesuffix='_agile_creation_spinup_' + f'{glacier_state}', delete=True)) print(f'{gdir.name} L mismatch at rgi_date:' f' {model.fls[-1].length_m - fls_spinup[-1].length_m} m') gdir.write_pickle(model.fls, 'model_flowlines', - filesuffix='_creation_spinup') + filesuffix='_creation_spinup_' + f'{glacier_state}') + + +def get_glacier_state_mb_model(gdir, glacier_state, save_to_gdir=False, + fls_spinup=None): + if fls_spinup is None: + fls_spinup = gdir.read_pickle('model_flowlines', + filesuffix='_creation_spinup_' + f'{glacier_state}') + model_type = None + model_args = {} + if glacier_state == 'equilibrium': + y0, halfsize = get_equilibrium_year_and_halfsize(gdir) + mb_model = MultipleFlowlineMassBalance( + gdir, + fls=fls_spinup, + mb_model_class=ConstantMassBalance, + y0=y0, + halfsize=halfsize, + filename='climate_historical', + ) + model_type = 'ConstantModel' + model_args['y0'] = y0 + model_args['halfsize'] = halfsize + model_args['filename'] = 'climate_historical' + + elif glacier_state == 'advancing': + t_bias_default = gdir.read_json('mb_calib')['temp_bias'] + repeat = False + ys = None + ye = None + if gdir.rgi_id == 'RGI60-16.02444': + t_bias_use = t_bias_default - 0.35 + repeat = True + ys = 1979 + ye = 2012 + elif gdir.rgi_id in ['RGI60-02.05098']: + t_bias_use = t_bias_default - 1.8 + elif gdir.rgi_id in ['RGI60-14.06794']: + t_bias_use = t_bias_default - 2.5 + elif gdir.rgi_id == 'RGI60-11.01450': + t_bias_use = t_bias_default - 2 + else: + raise NotImplementedError(f'{gdir.rgi_id}') + mb_model = MultipleFlowlineMassBalance( + gdir, + fls=fls_spinup, + mb_model_class=MonthlyTIModel, + temp_bias=t_bias_use, + repeat=repeat, + ys=ys, ye=ye, + filename='climate_historical') + + model_type = 'TIModel' + model_args['temp_bias'] = t_bias_use + model_args['repeat'] = repeat + model_args['ys'] = ys + model_args['ye'] = ye + model_args['filename'] = 'climate_historical' + + elif glacier_state == 'retreating': + t_bias_default = gdir.read_json('mb_calib')['temp_bias'] + if gdir.rgi_id == 'RGI60-16.02444': + t_bias_use = t_bias_default + 0 + elif gdir.rgi_id in ['RGI60-02.05098']: + t_bias_use = t_bias_default + 0 + elif gdir.rgi_id in ['RGI60-14.06794']: + t_bias_use = t_bias_default + 1.2 + elif gdir.rgi_id == 'RGI60-11.01450': + t_bias_use = t_bias_default + 1.3 + else: + raise NotImplementedError(f'{gdir.rgi_id}') + mb_model = MultipleFlowlineMassBalance( + gdir, + fls=fls_spinup, + mb_model_class=MonthlyTIModel, + temp_bias=t_bias_use, + filename='climate_historical') + + model_type = 'TIModel' + model_args['temp_bias'] = t_bias_use + model_args['filename'] = 'climate_historical' + + if save_to_gdir: + mb_models_agile = {'MB': {'type': model_type, + 'model_args': model_args, + 'years': np.array([1980, 2020])} + } + + # write the mb_model settings for later into gdir + gdir.write_pickle(mb_models_agile, 'inversion_input', + filesuffix=f'_agile_mb_models_{glacier_state}') + + return mb_model @entity_task(log, writes=['model_flowlines', 'inversion_input']) -def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, - yr_run_start, yr_run_end, +def evolve_glacier_and_create_measurements(gdir, glacier_state, used_mb_models, + yr_dmdt_start, yr_run_start, + yr_run_end, gcm='BCC-CSM2-MR', ssp='ssp370'): """TODO """ fls_spinup = gdir.read_pickle('model_flowlines', - filesuffix='_creation_spinup') + filesuffix='_creation_spinup_' + f'{glacier_state}') yr_rgi = gdir.rgi_date # use flow-parameters from default oggm inversion @@ -297,33 +473,26 @@ def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, # now start actual experiment run for the creation of measurements if used_mb_models == 'constant': - mb_models_agile = {'MB1': {'type': 'constant', - 'years': np.array([yr_run_start, - yr_dmdt_start])}, - 'MB2': {'type': 'constant', - 'years': np.array([yr_dmdt_start, - yr_run_end])} - } - mb_run = StackedMassBalance(gdir=gdir, - mb_model_settings=mb_models_agile) + ''' + mb_models_agile = {'MB1': {'type': 'constant', + 'years': np.array([yr_run_start, + yr_dmdt_start])}, + 'MB2': {'type': 'constant', + 'years': np.array([yr_dmdt_start, + yr_run_end])} + } + mb_run = StackedMassBalance(gdir=gdir, + mb_model_settings=mb_models_agile) + ''' + + raise NotImplementedError('Not implemented for different glacier' + 'states!') elif used_mb_models == 'TIModel': - mb_run = MultipleFlowlineMassBalance(gdir, - fls=fls_spinup, - mb_model_class=MonthlyTIModel, - filename='climate_historical') - - mb_models_agile = {'MB': {'type': 'TIModel', - 'years': []} - } - + mb_run = get_glacier_state_mb_model(gdir, glacier_state, + save_to_gdir=True) else: raise NotImplementedError(f'{used_mb_models}') - - # write the mb_model settings for later into gdir - gdir.write_pickle(mb_models_agile, 'inversion_input', - filesuffix='_agile_mb_models') - # do spinup period before first measurement model = SemiImplicitModel(copy.deepcopy(fls_spinup), mb_run, @@ -332,10 +501,12 @@ def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, model.run_until_and_store( yr_dmdt_start, diag_path=gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_dmdt_start', + filesuffix='_agile_true_dmdt_start_' + f'{glacier_state}', delete=True), fl_diag_path=gdir.get_filepath('fl_diagnostics', - filesuffix='_agile_true_dmdt_start', + filesuffix='_agile_true_dmdt_start_' + f'{glacier_state}', delete=True) ) @@ -348,14 +519,16 @@ def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, model.run_until_and_store( yr_rgi, diag_path=gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_init', + filesuffix='_agile_true_init_' + f'{glacier_state}', delete=True), fl_diag_path=gdir.get_filepath('fl_diagnostics', - filesuffix='_agile_true_init', + filesuffix='_agile_true_init_' + f'{glacier_state}', delete=True) ) gdir.write_pickle(model.fls, 'model_flowlines', - filesuffix='_agile_true_init') + filesuffix=f'_agile_true_init_{glacier_state}') rgi_date_area_km2 = model.area_km2 rgi_date_volume_km3 = model.volume_km3 rgi_date_fl_surface_h = model.fls[0].surface_h @@ -371,32 +544,35 @@ def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, model.run_until_and_store( yr_run_end, diag_path=gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_end', + filesuffix='_agile_true_end_' + f'{glacier_state}', delete=True), geom_path=gdir.get_filepath('model_geometry', - filesuffix='_agile_true_end', + filesuffix='_agile_true_end_' + f'{glacier_state}', delete=True), fl_diag_path=gdir.get_filepath('fl_diagnostics', - filesuffix='_agile_true_end', + filesuffix='_agile_true_end_' + f'{glacier_state}', delete=True) ) gdir.write_pickle(model.fls, 'model_flowlines', - filesuffix='_agile_true_end') + filesuffix=f'_agile_true_end_{glacier_state}') dmdtda_volume[1] = model.volume_m3 dmdtda_area[1] = model.area_m2 # agile model diagnostics to one file for the whole period tasks.merge_consecutive_run_outputs( gdir, - input_filesuffix_1='_agile_true_dmdt_start', - input_filesuffix_2='_agile_true_init', - output_filesuffix='_agile_true_total_run', + input_filesuffix_1=f'_agile_true_dmdt_start_{glacier_state}', + input_filesuffix_2=f'_agile_true_init_{glacier_state}', + output_filesuffix=f'_agile_true_total_run_{glacier_state}', delete_input=False) tasks.merge_consecutive_run_outputs( gdir, - input_filesuffix_1='_agile_true_total_run', - input_filesuffix_2='_agile_true_end', - output_filesuffix='_agile_true_total_run', + input_filesuffix_1=f'_agile_true_total_run_{glacier_state}', + input_filesuffix_2=f'_agile_true_end_{glacier_state}', + output_filesuffix=f'_agile_true_total_run_{glacier_state}', delete_input=False) # calculate dmdtda @@ -425,7 +601,7 @@ def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, 'fl_total_area:km2': {str(yr_rgi): rgi_date_area_km2}, } gdir.write_pickle(all_measurements, 'inversion_input', - filesuffix='_agile_measurements') + filesuffix=f'_agile_measurements_{glacier_state}') # conduct one future simulation for comparison # download locations for precipitation and temperature @@ -449,20 +625,22 @@ def evolve_glacier_and_create_measurements(gdir, used_mb_models, yr_dmdt_start, climate_input_filesuffix=rid, init_model_fls=model.fls, ys=yr_run_end, - output_filesuffix='_agile_true_future', + output_filesuffix='_agile_true_future_' + f'{glacier_state}', evolution_model=SemiImplicitModel, store_fl_diagnostics=True ) @entity_task(log, writes=['model_flowlines']) -def oggm_inversion_for_first_guess(gdir): +def oggm_inversion_for_first_guess(gdir, glacier_state): """TODO """ # now use _agile_true_init for an OGGM inversion for the first guess fls_inversion = gdir.read_pickle('inversion_flowlines')[0] fls_experiment = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] # load the oggm default dynamic parameters diag = gdir.get_diagnostics() @@ -474,7 +652,8 @@ def oggm_inversion_for_first_guess(gdir): ice_mask_index_diff = np.diff(np.where(fls_experiment.thick > 0.01)[0]) # if their is an ice-free gap raise an assertion assert np.all(ice_mask_index_diff == 1), (f'Their is a ice-free gap, check!' - f' (rgi_id: {gdir.rgi_id})') + f' (rgi_id: {gdir.rgi_id}, ' + f'glacier_state: {glacier_state})') # now use ice mask to create new inversion_flowlines fls_inversion.nx = int(sum(ice_mask)) @@ -542,19 +721,20 @@ def oggm_inversion_for_first_guess(gdir): assert np.allclose(fls_first_guess.widths_m, fls_experiment.widths_m) assert np.allclose(fls_first_guess.thick[ice_mask], inv_thick_use) gdir.write_pickle([fls_first_guess], 'model_flowlines', - filesuffix='_oggm_first_guess') + filesuffix=f'_oggm_first_guess_{glacier_state}') @entity_task(log, writes=['model_flowlines']) -def glabtop_inversion_for_first_guess(gdir): +def glabtop_inversion_for_first_guess(gdir, glacier_state): """TODO """ fls_experiment = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] ice_mask = np.where(fls_experiment.thick > 0.01, True, False) fg_thick = get_adaptive_upper_ice_thickness_limit( - fls_experiment, additional_ice_thickness=15) + fls_experiment, additional_ice_thickness=0) fg_bed_h = fls_experiment.surface_h[fls_experiment.thick > 0] - fg_thick bed_h_new = copy.deepcopy(fls_experiment.bed_h) @@ -579,11 +759,11 @@ def glabtop_inversion_for_first_guess(gdir): gdir=gdir) assert np.allclose(fl_fg.widths_m, fls_experiment.widths_m) - assert np.all(fl_fg._w0_m > 10) + assert np.all(fl_fg._w0_m >= 9.9) assert np.all(fl_fg.is_trapezoid) gdir.write_pickle([fl_fg], 'model_flowlines', - filesuffix='_glabtop_first_guess') + filesuffix=f'_glabtop_first_guess_{glacier_state}') def get_experiment_mb_model(gdir): @@ -594,27 +774,35 @@ def get_experiment_mb_model(gdir): return StackedMassBalance(gdir=gdir, mb_model_settings=mb_model_settings) -def oggm_default_initialisation(gdir, ys, ye, gcm='BCC-CSM2-MR', ssp='ssp370'): +def oggm_default_initialisation(gdir, glacier_state, ys, ye, + gcm='BCC-CSM2-MR', ssp='ssp370'): """ Use oggm default initialisation method and do a projection """ + # define mb_model depening on glacier_state + mb_model = get_glacier_state_mb_model(gdir, glacier_state) + # do dynamic initialisation spinup_model = tasks.run_dynamic_spinup( gdir, spinup_start_yr=ys, ye=ye, evolution_model=SemiImplicitModel, - model_flowline_filesuffix='_oggm_first_guess', - precision_absolute=0.1, - output_filesuffix='_oggm_dynamic_past', store_model_geometry=True, + model_flowline_filesuffix=f'_oggm_first_guess_{glacier_state}', + #precision_absolute=0.1, + mb_model_historical=mb_model, + output_filesuffix=f'_oggm_dynamic_past_{glacier_state}', + store_model_geometry=True, store_fl_diagnostics=True, store_model_evolution=True, ignore_errors=False, add_fixed_geometry_spinup=True) # do static initialisation fls_init = gdir.read_pickle('model_flowlines', - filesuffix='_oggm_first_guess') + filesuffix='_oggm_first_guess_' + f'{glacier_state}') static_model = tasks.run_from_climate_data( gdir, ye=ye, fixed_geometry_spinup_yr=ys, init_model_fls=fls_init, - output_filesuffix='_oggm_static_past', + mb_model=mb_model, + output_filesuffix=f'_oggm_static_past_{glacier_state}', store_model_geometry=True, store_fl_diagnostics=True) # conduct projection run for dynamic initialisation @@ -625,7 +813,8 @@ def oggm_default_initialisation(gdir, ys, ye, gcm='BCC-CSM2-MR', ssp='ssp370'): climate_input_filesuffix=rid, init_model_fls=spinup_model.fls, ys=ye, - output_filesuffix='_oggm_dynamic_future', + output_filesuffix='_oggm_dynamic_future_' + f'{glacier_state}', evolution_model=SemiImplicitModel, ) @@ -637,9 +826,10 @@ def oggm_default_initialisation(gdir, ys, ye, gcm='BCC-CSM2-MR', ssp='ssp370'): climate_input_filesuffix=rid, init_model_fls=static_model.fls, ys=ye, - output_filesuffix='_oggm_static_future', + output_filesuffix='_oggm_static_future_' + f'{glacier_state}', evolution_model=SemiImplicitModel, ) # calculate statistics - calculate_default_oggm_statistics(gdir) + calculate_default_oggm_statistics(gdir, glacier_state) diff --git a/agile1d/sandbox/define_idealized_experiment.py b/agile1d/sandbox/define_idealized_experiment.py index eadf522..b33b0c5 100644 --- a/agile1d/sandbox/define_idealized_experiment.py +++ b/agile1d/sandbox/define_idealized_experiment.py @@ -17,6 +17,7 @@ def idealized_experiment(use_experiment_glaciers, + use_experiment_glacier_states, inversion_settings_all, working_dir, output_folder, inversion_settings_individual=None, @@ -30,6 +31,13 @@ def idealized_experiment(use_experiment_glaciers, if override_params is None: override_params = {} + if any(state not in ['equilibrium', 'retreating', 'advancing'] + for state in use_experiment_glacier_states): + raise NotImplementedError("Glacier states could only be one of " + "['equilibrium', 'retreating', 'advancing']" + ", but given was " + f"{use_experiment_glacier_states}.") + utils.mkdir(working_dir) override_params['working_dir'] = working_dir @@ -56,6 +64,7 @@ def idealized_experiment(use_experiment_glaciers, 'L1-L2_files/elev_bands/' gdirs = create_idealized_experiments(use_experiment_glaciers, + use_experiment_glacier_states, prepro_border=cfg.PARAMS['border'], from_prepro_level=from_prepro_level, base_url=base_url, @@ -68,18 +77,27 @@ def idealized_experiment(use_experiment_glaciers, all_experiments = [] for inv_setting in inversion_settings_all: for gdir in gdirs: - inv_setting_use = copy.deepcopy(inv_setting) - if inversion_settings_individual is not None: - try: - individual_setting = inversion_settings_individual[gdir.rgi_id] - for key in individual_setting.keys(): - inv_setting_use[key] = individual_setting[key] - except KeyError: - # if no individual settings for this glacier move on - pass - all_experiments.append((gdir, dict(inversion_settings=inv_setting_use, - output_folder=output_folder) - )) + for glacier_state in use_experiment_glacier_states: + inv_setting_use = copy.deepcopy(inv_setting) + if inversion_settings_individual is not None: + try: + individual_setting = \ + inversion_settings_individual[gdir.rgi_id] + for key in individual_setting.keys(): + inv_setting_use[key] = individual_setting[key] + except KeyError: + # if no individual settings for this glacier move on + pass + # add glacier state to outputfilesuffix + inv_setting_use['experiment_description'] = ( + f"{glacier_state}_" + f"{inv_setting_use['experiment_description']}") + all_experiments.append((gdir, + dict( + glacier_state=glacier_state, + inversion_settings=inv_setting_use, + output_folder=output_folder) + )) workflow.execute_entity_task(conduct_sandbox_inversion, all_experiments, init_model_fls=init_model_fls, @@ -116,23 +134,39 @@ def add_future_projection_run(gdir, data_logger, gcm='BCC-CSM2-MR', @entity_task(log, writes=['inversion_input', 'model_flowlines']) -def conduct_sandbox_inversion(gdir, inversion_settings=None, +def conduct_sandbox_inversion(gdir, glacier_state, inversion_settings=None, output_folder=None, init_model_fls='_oggm_first_guess', gcm='BCC-CSM2-MR', ssp='ssp370', - print_statistic=True): + print_statistic=True, + return_data_logger=False): """TODO""" - # check if mb_model_settings should be loaded from gdir - if inversion_settings['mb_models_settings'] == 'load_from_gdir': - inversion_settings['mb_models_settings'] = \ - gdir.read_pickle('inversion_input', - filesuffix='_agile_mb_models') + # add glacier state to first guess + init_model_fls = f'{init_model_fls}_{glacier_state}' + + # load mb_model settings from gdir depending on glacier state + inversion_settings['mb_models_settings'] = \ + gdir.read_pickle('inversion_input', + filesuffix='_agile_mb_models_' + f'{glacier_state}') + + # if perfect spinup we must add glacier state to string + spn_option = list(inversion_settings['spinup_options'].keys())[0] + if spn_option in ['perfect_sfc_h', 'perfect_thickness', 'perfect_section']: + inversion_settings['spinup_options'][spn_option] = \ + inversion_settings['spinup_options'][spn_option] + f'_{glacier_state}' + + # also for perfect_bed_h we need to add the glacier state + if 'perfect_bed_h' in list(inversion_settings['spinup_options'].keys()): + inversion_settings['spinup_options']['perfect_bed_h'] = \ + inversion_settings['spinup_options']['perfect_bed_h'] + f'_{glacier_state}' # load the observations to use from gdir all_measurements = gdir.read_pickle('inversion_input', - filesuffix='_agile_measurements') + filesuffix='_agile_measurements_' + f'{glacier_state}') used_measurements = copy.deepcopy(inversion_settings['observations']) inversion_settings['observations'] = {} for msr in used_measurements: @@ -144,20 +178,21 @@ def conduct_sandbox_inversion(gdir, inversion_settings=None, else: raise RuntimeError(f'No observation for {msr} available!') - prepare_for_agile_inversion(gdir, - inversion_settings=inversion_settings, - filesuffix= - inversion_settings['experiment_description']) + inversion_filesuffix = f"_{inversion_settings['experiment_description']}" + + prepare_for_agile_inversion( + gdir, inversion_settings=inversion_settings, + filesuffix=inversion_filesuffix) data_logger = agile_inversion(gdir, inversion_input_filesuffix= - inversion_settings['experiment_description'], + inversion_filesuffix, init_model_filesuffix=None, init_model_fls=init_model_fls, climate_filename='climate_historical', climate_filesuffix='', - output_filesuffix='_agile_output_' + - inversion_settings['experiment_description'], + output_filesuffix= + f'_agile_output{inversion_filesuffix}', output_filepath=output_folder, save_dataset=True, give_data_logger_back=True) @@ -165,5 +200,9 @@ def conduct_sandbox_inversion(gdir, inversion_settings=None, add_future_projection_run(gdir, data_logger=data_logger, gcm=gcm, ssp=ssp) calculate_result_statistics(gdir, + glacier_state=glacier_state, data_logger=data_logger, print_statistic=print_statistic) + + if return_data_logger: + return data_logger diff --git a/agile1d/sandbox/graphics.py b/agile1d/sandbox/graphics.py index 732682c..8c13538 100644 --- a/agile1d/sandbox/graphics.py +++ b/agile1d/sandbox/graphics.py @@ -1,5 +1,6 @@ from oggm import cfg, utils -from agile1d.sandbox.create_glaciers_with_measurements import create_idealized_experiments +from agile1d.sandbox.create_glaciers_with_measurements import create_idealized_experiments, \ + get_equilibrium_year_and_halfsize, get_spinup_mb_model, get_glacier_state_mb_model import matplotlib.pyplot as plt from matplotlib import colors as m_colors import numpy as np @@ -10,7 +11,8 @@ import os from oggm import graphics from oggm.graphics import _plot_map, truncate_colormap -from oggm.core.massbalance import MonthlyTIModel, MultipleFlowlineMassBalance, RandomMassBalance +from oggm.core.massbalance import MonthlyTIModel, MultipleFlowlineMassBalance, RandomMassBalance, \ + ConstantMassBalance from oggm.core.flowline import FlowlineModel import salem from oggm.core import gis @@ -134,14 +136,15 @@ def plot_outline_with_dem(gdir, ax=None, map_extend_factor=3, use_netcdf=False): # heading text -def plot_heading_text(gdir, fig, ax, parameters_text_height=0.5, linespacing=1): +def plot_heading_text(gdir, glacier_state, fig, ax, parameters_text_height=0.5, + linespacing=1): ax.axis('off') ax.text( 0.5, 1, - f'{gdir.rgi_id}: {gdir.name}', + f'{gdir.rgi_id}: {gdir.name}, {glacier_state}', fontsize=16, horizontalalignment='center', verticalalignment='top', @@ -174,16 +177,24 @@ def plot_heading_text(gdir, fig, ax, parameters_text_height=0.5, linespacing=1): # glacier bed with hypsomentry # adapted from graphics.plot_modeloutput_section -def plot_bed_hypsometry(gdir, ax, fig, grid_points_added=5, plot_non_trap=True, - bins=20, legend_x_position=0.8, legend_y_position=1.0, - add_future=True, ye=2050): +def plot_bed_hypsometry(gdir, glacier_state, ax, fig, grid_points_added=5, + plot_non_trap=True, bins=20, legend_x_position=0.8, + legend_y_position=1.0, add_future=True, ye=2050): # open flowlines to be plotted (only for single flowline models) - fls_true_start = gdir.read_pickle('model_flowlines', filesuffix='_creation_spinup')[0] - fls_true_rgi = gdir.read_pickle('model_flowlines', filesuffix='_agile_true_init')[0] - fls_true_end = gdir.read_pickle('model_flowlines', filesuffix='_agile_true_end')[0] + fls_true_start = gdir.read_pickle('model_flowlines', + filesuffix='_creation_spinup_' + f'{glacier_state}')[0] + fls_true_rgi = gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] + fls_true_end = gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_end_' + f'{glacier_state}')[0] if add_future: - fls_true_future = get_fl_diagnostics(gdir, '_agile_true_future') + fls_true_future = get_fl_diagnostics(gdir, + '_agile_true_future_' + f'{glacier_state}') fls_true_future = fls_true_future.loc[{'time': ye}] fls_true_future_thickness = fls_true_future.thickness_m.values fls_true_future_surface_h = fls_true_future.bed_h.values + fls_true_future_thickness @@ -315,7 +326,7 @@ def plot_glacier_surface(fl, label, color): # Mass Balance plot -def plot_specific_mb(gdir, ax, +def plot_specific_mb(gdir, glacier_state, ax, ys=1920, ye=2100, y_switch_spinup_run=1980, @@ -325,17 +336,11 @@ def plot_specific_mb(gdir, ax, # define mass balance models to use fls_spinup = gdir.read_pickle('model_flowlines', filesuffix='_consensus') - mb_model_random = MultipleFlowlineMassBalance(gdir, - fls=fls_spinup, - mb_model_class=RandomMassBalance, - seed=2, - filename='climate_historical', - y0=1950, - halfsize=30) - mb_model_run = MultipleFlowlineMassBalance(gdir, - fls=fls_spinup, - mb_model_class=MonthlyTIModel, - filename='climate_historical') + + mb_model_random = get_spinup_mb_model(gdir, glacier_state, fls_spinup) + mb_model_run = get_glacier_state_mb_model(gdir, glacier_state, + fls_spinup=fls_spinup) + gcm = 'BCC-CSM2-MR' ssp = 'ssp370' rid = '_{}_{}'.format(gcm, ssp) @@ -345,7 +350,8 @@ def plot_specific_mb(gdir, ax, filename='gcm_data', input_filesuffix=rid) fls = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init') + filesuffix='_agile_true_init_' + f'{glacier_state}') # read out mb values specific_mb = [] @@ -400,34 +406,39 @@ def plot_specific_mb(gdir, ax, # glacier evolution during creation -def plot_glacier_evolution(gdir, ax_v, ax_a, ax_l, text_position=1.1, - add_text=True, add_future=True, add_fg=True, - ye=2100): +def plot_glacier_evolution(gdir, glacier_state, ax_v, ax_a, ax_l, + text_position=1.1, add_text=True, add_future=True, + add_fg=True, ye=2100): color = 'C0' color_fg = 'C1' # open glacier creation with xr.open_dataset(gdir.get_filepath('model_diagnostics', - filesuffix='_agile_creation_spinup')) as ds: + filesuffix='_agile_creation_spinup_' + f'{glacier_state}')) as ds: ds_spinup = ds.load() # open actual historical run with xr.open_dataset(gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_total_run')) as ds: + filesuffix='_agile_true_total_run_' + f'{glacier_state}')) as ds: ds_run = ds.load() # open future run with xr.open_dataset(gdir.get_filepath('model_diagnostics', - filesuffix='_agile_true_future')) as ds: + filesuffix='_agile_true_future_' + f'{glacier_state}')) as ds: ds_future = ds.load() if add_fg: with xr.open_dataset(gdir.get_filepath('model_diagnostics', - filesuffix='_oggm_dynamic_past')) as ds: + filesuffix='_oggm_static_past_' + f'{glacier_state}')) as ds: ds_run_fg = ds.load() with xr.open_dataset(gdir.get_filepath('model_diagnostics', - filesuffix='_oggm_dynamic_future')) as ds: + filesuffix='_oggm_static_future_' + f'{glacier_state}')) as ds: ds_future_fg = ds.load() if add_future: @@ -500,10 +511,15 @@ def plot_glacier_evolution(gdir, ax_v, ax_a, ax_l, text_position=1.1, # dhdt, bed diff plot -def plot_dhdt_first_guess_db(gdir, ax_dh, ax_db, ax_ds=None, grid_points_added=0): +def plot_dhdt_first_guess_db(gdir, glacier_state, ax_dh, ax_db, ax_ds=None, + grid_points_added=0): # open flowlines to be plotted (only for single flowline models) - fls_true_rgi = gdir.read_pickle('model_flowlines', filesuffix='_agile_true_init')[0] - fls_first_guess = gdir.read_pickle('model_flowlines', filesuffix='_oggm_first_guess')[0] + fls_true_rgi = gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] + fls_first_guess = gdir.read_pickle('model_flowlines', + filesuffix='_oggm_first_guess_' + f'{glacier_state}')[0] # define extend to plot according to longest flowline @@ -514,10 +530,14 @@ def plot_dhdt_first_guess_db(gdir, ax_dh, ax_db, ax_ds=None, grid_points_added=0 max_grid_point += grid_points_added # open different flowline diagnostics - ds_start = get_fl_diagnostics(gdir, '_agile_true_dmdt_start') - ds_init = get_fl_diagnostics(gdir, '_agile_true_init') - ds_end = get_fl_diagnostics(gdir, '_agile_true_end') - ds_default = get_fl_diagnostics(gdir, '_oggm_dynamic_past').sel({'time': 2000}) + ds_start = get_fl_diagnostics(gdir, '_agile_true_dmdt_start_' + f'{glacier_state}') + ds_init = get_fl_diagnostics(gdir, '_agile_true_init_' + f'{glacier_state}') + ds_end = get_fl_diagnostics(gdir, '_agile_true_end_' + f'{glacier_state}') + ds_default = get_fl_diagnostics(gdir, '_oggm_static_past_' + f'{glacier_state}').sel({'time': 2000}) ds_all = xr.merge([ds_start, ds_init, ds_end]) @@ -580,11 +600,17 @@ def plot_dhdt_first_guess_db(gdir, ax_dh, ax_db, ax_ds=None, grid_points_added=0 # plot velocities -def plot_velocities(gdir, ax, grid_points_added=5): +def plot_velocities(gdir, glacier_state, ax, grid_points_added=5): # open flowlines to be plotted (only for single flowline models) - fls_true_start = gdir.read_pickle('model_flowlines', filesuffix='_creation_spinup')[0] - fls_true_rgi = gdir.read_pickle('model_flowlines', filesuffix='_agile_true_init')[0] - fls_true_end = gdir.read_pickle('model_flowlines', filesuffix='_agile_true_end')[0] + fls_true_start = gdir.read_pickle('model_flowlines', + filesuffix='_creation_spinup_' + f'{glacier_state}')[0] + fls_true_rgi = gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] + fls_true_end = gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_end_' + f'{glacier_state}')[0] # define extend to plot according to longest flowline max_grid_point = np.max( @@ -594,9 +620,12 @@ def plot_velocities(gdir, ax, grid_points_added=5): max_grid_point += grid_points_added # open different flowline diagnostics - ds_start = get_fl_diagnostics(gdir, '_agile_true_dmdt_start') - ds_init = get_fl_diagnostics(gdir, '_agile_true_init') - ds_end = get_fl_diagnostics(gdir, '_agile_true_end') + ds_start = get_fl_diagnostics(gdir, '_agile_true_dmdt_start_' + f'{glacier_state}') + ds_init = get_fl_diagnostics(gdir, '_agile_true_init_' + f'{glacier_state}') + ds_end = get_fl_diagnostics(gdir, '_agile_true_end_' + f'{glacier_state}') ds_all = xr.merge([ds_start, ds_init, ds_end]) @@ -620,8 +649,8 @@ def plot_velocities(gdir, ax, grid_points_added=5): # create the whole plot -def create_synthetic_glacier_plot(gdir, fig, save_fig=False, fig_dir='', - legend_position=(0.2, 1.2), +def create_synthetic_glacier_plot(gdir, glacier_state, fig, save_fig=False, + fig_dir='', legend_position=(0.2, 1.2), map_extend_factor=2, add_future=False, add_fg=True, ye=2020): top_margin = 0.02 @@ -700,16 +729,17 @@ def create_synthetic_glacier_plot(gdir, fig, save_fig=False, fig_dir='', # heading plots plot_outline_with_dem(gdir, ax=ax_map, map_extend_factor=map_extend_factor) - plot_heading_text(gdir, fig, ax_text, + plot_heading_text(gdir, glacier_state, fig, ax_text, parameters_text_height=0.7, linespacing=1.3) # left column plots - plot_bed_hypsometry(gdir, ax=ax_bed, fig=fig, grid_points_added=15, bins=15, + plot_bed_hypsometry(gdir, glacier_state, ax=ax_bed, fig=fig, + grid_points_added=15, bins=15, legend_x_position=legend_position[0], legend_y_position=legend_position[1], add_future=add_future, ye=ye) - plot_dhdt_first_guess_db(gdir, + plot_dhdt_first_guess_db(gdir, glacier_state, ax_dh=ax_dh, ax_db=ax_db, # ax_ds=ax_ds, @@ -730,9 +760,10 @@ def create_synthetic_glacier_plot(gdir, fig, save_fig=False, fig_dir='', ax_dh.set_xticklabels(labels) # right column plots - plot_specific_mb(gdir, ax=ax_mb, ys=1920, ye=ye, + plot_specific_mb(gdir, glacier_state, ax=ax_mb, ys=1920, ye=ye, add_text=True, add_future=add_future) plot_glacier_evolution(gdir, + glacier_state, ax_v=ax_v, ax_a=ax_a, ax_l=None, @@ -768,7 +799,15 @@ def create_synthetic_glacier_plot(gdir, fig, save_fig=False, fig_dir='', fig.savefig(os.path.join(result_folder, f'{gdir.name}.png')) -def create_and_save_all_synthetic_glacier_plots(fig_dir='', **kwargs): +def create_and_save_all_synthetic_glacier_plots(fig_dir='', + glacier_names=('Baltoro', + 'Aletsch', + 'Artesonraju', + 'Peyto'), + glacier_states=('equilibrium', + 'retreating', + 'advancing'), + **kwargs): if fig_dir == '': raise ValueError('You should define a directory where to save the ' @@ -784,14 +823,9 @@ def create_and_save_all_synthetic_glacier_plots(fig_dir='', **kwargs): base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/' \ 'L1-L2_files/elev_bands/' - glacier_names = ['Baltoro', - 'Aletsch', - 'Artesonraju', - 'Peyto', - ] - cfg.PARAMS['cfl_number'] = 0.5 gdirs = create_idealized_experiments(glacier_names, + glacier_states, prepro_border=prepro_border, from_prepro_level=from_prepro_level, base_url=base_url, ) @@ -811,14 +845,16 @@ def create_and_save_all_synthetic_glacier_plots(fig_dir='', **kwargs): } for gdir in gdirs: - kwargs_use = glacier_kwargs[gdir.rgi_id] - save_fig = True + for glacier_state in glacier_states: + kwargs_use = glacier_kwargs[gdir.rgi_id] + save_fig = True - fig = plt.figure(figsize=(10, 6)) + fig = plt.figure(figsize=(10, 6)) - for key, item in kwargs.items(): - kwargs_use[key] = item + for key, item in kwargs.items(): + kwargs_use[key] = item - create_synthetic_glacier_plot(gdir, fig, save_fig=save_fig, - fig_dir=fig_dir, - **kwargs_use) + create_synthetic_glacier_plot(gdir, glacier_state, fig, + save_fig=save_fig, + fig_dir=fig_dir, + **kwargs_use) diff --git a/agile1d/sandbox/individual_experiment_dashboard.py b/agile1d/sandbox/individual_experiment_dashboard.py index 168a2d9..824c929 100644 --- a/agile1d/sandbox/individual_experiment_dashboard.py +++ b/agile1d/sandbox/individual_experiment_dashboard.py @@ -33,6 +33,7 @@ def find_class(self, module, name): all_experiment_settings = {} experiment_select = [] glacier_select = [] +glacier_state_select = [] menu = pn.Column() button = pn.widgets.Button(name='Select new one', button_type='primary') @@ -59,6 +60,7 @@ def define_options_for_experiment(event): options_selects.append(pn.widgets.Select(name=name_tmp, options=options_tmp)) new_menu = pn.Column(glacier_select, + glacier_state_select, experiment_select) for opt_select in options_selects: @@ -75,6 +77,7 @@ def individual_experiment_dashboard(working_dir, input_folder, global all_experiment_settings global experiment_select global glacier_select + global glacier_state_select cfg.PATHS['working_dir'] = working_dir gdirs = workflow.init_glacier_directories() @@ -82,6 +85,7 @@ def individual_experiment_dashboard(working_dir, input_folder, # Extract info from directory all_files = os.listdir(input_folder) all_glaciers = [] + all_glacier_states = [] all_experiments = [] # only keep experiment files @@ -99,11 +103,16 @@ def individual_experiment_dashboard(working_dir, input_folder, if glacier not in all_glaciers: all_glaciers.append(glacier) + # second is always glacier state + single_glacier_state = splits[1] + if single_glacier_state not in all_glacier_states: + all_glacier_states.append(single_glacier_state) + # loop through the rest, if three characters followed by numbers it is # an experiment setting experiment = '' experiment_settings_tmp = [] - for single_split in splits[1:]: + for single_split in splits[2:]: # find out if it has 3 characters followed only by numbers, if so # is is an experiment setting if len(single_split) > 3: @@ -144,6 +153,8 @@ def AERR(a1, a2): # Define menu glacier_select = pn.widgets.Select(name='glacier', options=all_glaciers) + glacier_state_select = pn.widgets.Select(name='glacier_state', + options=all_glacier_states) experiment_select = pn.widgets.Select(name='experiment', options=all_experiments) experiment_select.param.watch(define_options_for_experiment, 'value') @@ -163,7 +174,7 @@ def get_description_accordion(): return accordion # Define individual plots - def get_individual_plot(open_file, current_file): + def get_individual_plot(open_file, current_file, glacier_state): ds = open_file # get reference flowline for true values @@ -172,7 +183,8 @@ def get_individual_plot(open_file, current_file): for gdir in gdirs: if gdir.rgi_id == rgi_id: fl_ref = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] # now calculate data for delta bed_h and w0_m data_bed_h = [] @@ -426,11 +438,14 @@ def get_heatmap(data, lim, title, kdim='x', vdim='Iteration', height=200): for gdir in gdirs: if gdir.rgi_id == rgi_id: fl_ref_rgi = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] fl_ref_start = gdir.read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] + filesuffix='_creation_spinup_' + f'{glacier_state}')[0] fl_ref_end = gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_end')[0] + filesuffix='_agile_true_end_' + f'{glacier_state}')[0] def get_performance_sfc_h_array(fct, data, ref_val): return [np.around(fct(val, ref_val), decimals=2) for val in data] @@ -770,6 +785,7 @@ def get_curve(x, y, label, color, height=200): # put together and link menu with plots current_file_first = list(compress(all_files, [glacier_select.value in file and + glacier_state_select.value in file and experiment_select.value in file for file in all_files])) for opt_select in options_selects: @@ -784,12 +800,14 @@ def get_curve(x, y, label, color, height=200): open_file = CpuUnpickler(handle).load() # pickle.load(handle) - figure = get_individual_plot(open_file, current_file_first) + figure = get_individual_plot(open_file, current_file_first, + glacier_state_select.value) def change_figure(event): # here get the right filename for the current selection current_file = list(compress(all_files, [glacier_select.value + '_' + + glacier_state_select.value + '_' + experiment_select.value in file for file in all_files])) @@ -816,7 +834,8 @@ def change_figure(event): open_file = CpuUnpickler(handle).load() # pickle.load(handle,) - figure.objects = [get_individual_plot(open_file, current_file)] + figure.objects = [get_individual_plot(open_file, current_file, + glacier_state_select.value)] button.on_click(change_figure) diff --git a/agile1d/sandbox/run_idealized_experiment.py b/agile1d/sandbox/run_idealized_experiment.py index 44584bc..05ba2fb 100644 --- a/agile1d/sandbox/run_idealized_experiment.py +++ b/agile1d/sandbox/run_idealized_experiment.py @@ -70,10 +70,14 @@ def parse_args(args): if not hasattr(foo, 'use_experiment_glaciers'): raise InvalidParamsError('No variable "experiment_glaciers" in ' 'provided experiments file!') + if not hasattr(foo, 'use_experiment_glacier_states'): + raise InvalidParamsError('No variable "experiment_glacier_states" in ' + 'provided experiments file!') if not hasattr(foo, 'inversion_settings_individual'): raise InvalidParamsError('No variable "inversion_settings_individual" in ' 'provided experiments file!') use_experiment_glaciers = foo.use_experiment_glaciers + use_experiment_glacier_states = foo.use_experiment_glacier_states for glacier in use_experiment_glaciers: if glacier not in experiment_glaciers.keys(): raise InvalidParamsError(f'{glacier} not supported! First must be ' @@ -86,6 +90,7 @@ def parse_args(args): inversion_settings_individual=foo.inversion_settings_individual, init_model_fls=init_model_fls, use_experiment_glaciers=use_experiment_glaciers, + use_experiment_glacier_states=use_experiment_glacier_states, print_statistic=args.print_statistic ) diff --git a/agile1d/tests/test_sandbox.py b/agile1d/tests/test_sandbox.py index 463ecbb..30e8579 100644 --- a/agile1d/tests/test_sandbox.py +++ b/agile1d/tests/test_sandbox.py @@ -6,21 +6,25 @@ import os import pickle +from agile1d.core.cost_function import initialise_mb_models from agile1d.core.inversion import get_default_inversion_settings from agile1d.core.massbalance import StackedMassBalance -from agile1d.sandbox.create_glaciers_with_measurements import create_idealized_experiments -from agile1d.sandbox.define_idealized_experiment import idealized_experiment +from agile1d.sandbox.create_glaciers_with_measurements import create_idealized_experiments, \ + get_glacier_state_mb_model +from agile1d.sandbox.define_idealized_experiment import idealized_experiment, \ + conduct_sandbox_inversion from oggm import cfg +from oggm.core.massbalance import ConstantMassBalance, MonthlyTIModel do_plot = False pytestmark = [pytest.mark.filterwarnings("ignore: " - "should not be instantiated.:DeprecationWarning"), + "should not be instantiated.:DeprecationWarning"), pytest.mark.test_env("sandbox")] class TestSandbox: - def test_create_idealized_experiments(self, test_dir): + def test_create_idealized_experiments_and_mb_glacier_state(self, test_dir): cfg.initialize(logging_level='WARNING') cfg.PATHS['working_dir'] = test_dir cfg.PARAMS['use_multiprocessing'] = False @@ -39,10 +43,12 @@ def test_create_idealized_experiments(self, test_dir): 'Artesonraju', 'Peyto' ] + glacier_states = ['equilibrium', 'retreating', 'advancing'] cfg.PARAMS['use_multiprocessing'] = False cfg.PARAMS['cfl_number'] = 0.5 gdirs = create_idealized_experiments(glacier_names, + glacier_states=glacier_states, prepro_border=prepro_border, from_prepro_level=from_prepro_level, base_url=base_url, ) @@ -63,41 +69,54 @@ def resulting_run_test(filesuffix, ys, ye): # test that the file of the creation spinup exist for gdir in gdirs: - gdir.has_file('model_diagnostics', - filesuffix='_agile_creation_spinup') + for glacier_state in glacier_states: + gdir.has_file('model_diagnostics', + filesuffix='_agile_creation_spinup_' + f'{glacier_state}') # test that resulting true run dataset contain the complete period for gdir in gdirs: - resulting_run_test(filesuffix='_agile_true_total_run', - ys=1980, - ye=2020) + for glacier_state in glacier_states: + resulting_run_test(filesuffix='_agile_true_total_run_' + f'{glacier_state}', + ys=1980, + ye=2020) # test that the resulting gdirs contain the future climate file for gdir in gdirs: assert gdir.has_file('gcm_data', filesuffix='_BCC-CSM2-MR_ssp370') - resulting_run_test(filesuffix='_agile_true_future', - ys=2020, - ye=2101) + for glacier_state in glacier_states: + resulting_run_test(filesuffix='_agile_true_future_' + f'{glacier_state}', + ys=2020, + ye=2101) # test that the resulting gdirs contain the oggm default run files for gdir in gdirs: - resulting_run_test(filesuffix='_oggm_dynamic_past', - ys=1980, - ye=2020) - assert gdir.has_file('fl_diagnostics', - filesuffix='_oggm_dynamic_past') - resulting_run_test(filesuffix='_oggm_dynamic_future', - ys=2020, - ye=2101) - - resulting_run_test(filesuffix='_oggm_static_past', - ys=1980, - ye=2020) - assert gdir.has_file('fl_diagnostics', - filesuffix='_oggm_static_past') - resulting_run_test(filesuffix='_oggm_static_future', - ys=2020, - ye=2101) + for glacier_state in glacier_states: + resulting_run_test(filesuffix='_oggm_dynamic_past_' + f'{glacier_state}', + ys=1980, + ye=2020) + assert gdir.has_file('fl_diagnostics', + filesuffix='_oggm_dynamic_past_' + f'{glacier_state}') + resulting_run_test(filesuffix='_oggm_dynamic_future_' + f'{glacier_state}', + ys=2020, + ye=2101) + + resulting_run_test(filesuffix='_oggm_static_past_' + f'{glacier_state}', + ys=1980, + ye=2020) + assert gdir.has_file('fl_diagnostics', + filesuffix='_oggm_static_past_' + f'{glacier_state}') + resulting_run_test(filesuffix='_oggm_static_future_' + f'{glacier_state}', + ys=2020, + ye=2101) # test that the resulting gdirs contain the oggm default statistics def all_stats_finite(ds, use_year, is_static=False): @@ -116,113 +135,176 @@ def all_stats_finite(ds, use_year, is_static=False): assert np.all(np.isfinite(ds[var][year])) for gdir in gdirs: - fp = os.path.join(gdir.dir, - 'default_oggm_statistics.pkl') - with open(fp, 'rb') as handle: - ds_default_stats = pickle.load(handle) - - for stat_key in ds_default_stats.keys(): - stat_suffixes = stat_key.split('_')[-1] - pure_key = stat_key.removesuffix('_' + stat_suffixes) - if pure_key in ['observations_stats', 'controls_stats', - 'past_evol_stats', 'today_state_stats', - 'future_evol_stats', 'past_state_stats']: - use_year = False - if pure_key in ['observations_stats']: - use_year = True - all_stats_finite(ds_default_stats[stat_key], use_year, - stat_suffixes=='static') - else: - raise NotImplementedError(f'{stat_key}') + for glacier_state in glacier_states: + fp = os.path.join(gdir.dir, + f'default_oggm_statistics_{glacier_state}.pkl') + with open(fp, 'rb') as handle: + ds_default_stats = pickle.load(handle) + + for stat_key in ds_default_stats.keys(): + stat_suffixes = stat_key.split('_')[-1] + pure_key = stat_key.removesuffix('_' + stat_suffixes) + if pure_key in ['observations_stats', 'controls_stats', + 'past_evol_stats', 'today_state_stats', + 'future_evol_stats', 'past_state_stats']: + use_year = False + if pure_key in ['observations_stats']: + use_year = True + all_stats_finite(ds_default_stats[stat_key], use_year, + stat_suffixes == 'static') + else: + raise NotImplementedError(f'{stat_key}') - assert gdir.has_file('model_flowlines', - filesuffix='_oggm_first_guess') - assert gdir.has_file('model_flowlines', - filesuffix='_glabtop_first_guess') + assert gdir.has_file('model_flowlines', + filesuffix='_oggm_first_guess_' + f'{glacier_state}') + assert gdir.has_file('model_flowlines', + filesuffix='_glabtop_first_guess_' + f'{glacier_state}') if do_plot: for gdir in gdirs: - fl_oggm = gdir.read_pickle('model_flowlines')[0] - fl_consensus = gdir.read_pickle('model_flowlines', - filesuffix='_consensus')[0] - fl_spinup = gdir.read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] - fl_agile_init = \ - gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_init')[0] - fl_agile_end = \ - gdir.read_pickle('model_flowlines', - filesuffix='_agile_true_end')[0] - fl_oggm_first_guess = \ - gdir.read_pickle('model_flowlines', - filesuffix='_oggm_first_guess')[0] - - def get_fl_diagnostics(filesuffix): - f = gdir.get_filepath('fl_diagnostics', - filesuffix=filesuffix) - with xr.open_dataset(f, group=f'fl_0') as ds: - ds = ds.load() - return ds - - fl_diag_init = get_fl_diagnostics('_agile_true_dmdt_start') - fl_diag_start = get_fl_diagnostics('_agile_true_init') - fl_diag_end = get_fl_diagnostics('_agile_true_end') - - fig = plt.figure(figsize=(15, 10)) - (ax, ax2, ax3, ax4, ax5, ax6) = fig.subplots(6, 1) - # ax.plot(fl_oggm.dis_on_line, fl_oggm.bed_h, - # label='OGGM Flowline') - ax.plot(fl_spinup.dis_on_line, fl_spinup.bed_h, - label='Spinup Flowline') - ax.plot(fl_oggm.dis_on_line, fl_oggm.surface_h, - label='OGGM surface_h') - ax.plot(fl_spinup.dis_on_line, fl_spinup.surface_h, - label='Spinup surface_h') - # ax.plot(fl_consensus.dis_on_line, fl_consensus.surface_h, - # label='Consensus surface_h') - ax.plot(fl_agile_init.dis_on_line, fl_agile_init.surface_h, - label='Init surface_h') - ax.plot(fl_agile_end.dis_on_line, fl_agile_end.surface_h, - label='END surface_h') - ax.legend(); - - ax2.axhline(0, color='black') - ax2.plot(fl_spinup.surface_h - fl_agile_init.surface_h, - label='positive means retreating ' - '(spinup to rgi _date)') - ax2.plot(fl_agile_init.surface_h - fl_agile_end.surface_h, - label='rgi_date to end') - ax2.set_title('delta surface_h') - ax2.legend() - - ax3.axhline(0, color='black') - ax3.plot(fl_spinup.bed_h - fl_oggm_first_guess.bed_h, - label='positive means overdeepening') - ax3.set_title('delta bed_h') - ax3.legend() - - ax4.plot(fl_oggm_first_guess.is_trapezoid, label='trapez') - ax4.plot(fl_oggm_first_guess.is_rectangular, label='rect') - ax4.legend() - - ax5.plot(fl_oggm.thick > 0, label='OGGM thick > 0') - ax5.plot(fl_agile_init.thick > 0, - label='agile init thick > 0') - ax5.legend() - - ax6.plot(fl_diag_init.ice_velocity_myr[-1], - label=f'year {fl_diag_init.time[-1].values}') - ax6.plot(fl_diag_start.ice_velocity_myr[-1], - label=f'year {fl_diag_start.time[-1].values} ' - f'(rgi date)') - ax6.plot(fl_diag_end.ice_velocity_myr[-1], - label=f'year {fl_diag_end.time[-1].values}') - ax6.set_title('velocity') - ax6.legend() - - fig.suptitle(gdir.name, fontsize=16) - fig.tight_layout(pad=1.0) - plt.show() + for glacier_state in glacier_states: + fl_oggm = gdir.read_pickle('model_flowlines')[0] + fl_consensus = gdir.read_pickle('model_flowlines', + filesuffix='_consensus')[0] + fl_spinup = gdir.read_pickle('model_flowlines', + filesuffix='_creation_spinup' + f'_{glacier_state}')[0] + fl_agile_init = \ + gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_init_' + f'{glacier_state}')[0] + fl_agile_end = \ + gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_end_' + f'{glacier_state}')[0] + fl_oggm_first_guess = \ + gdir.read_pickle('model_flowlines', + filesuffix='_oggm_first_guess_' + f'{glacier_state}')[0] + + def get_fl_diagnostics(filesuffix): + f = gdir.get_filepath('fl_diagnostics', + filesuffix=filesuffix) + with xr.open_dataset(f, group=f'fl_0') as ds: + ds = ds.load() + return ds + + fl_diag_init = get_fl_diagnostics('_agile_true_dmdt_start_' + f'{glacier_state}') + fl_diag_start = get_fl_diagnostics('_agile_true_init_' + f'{glacier_state}') + fl_diag_end = get_fl_diagnostics('_agile_true_end_' + f'{glacier_state}') + + fig = plt.figure(figsize=(15, 10)) + (ax, ax2, ax3, ax4, ax5, ax6) = fig.subplots(6, 1) + # ax.plot(fl_oggm.dis_on_line, fl_oggm.bed_h, + # label='OGGM Flowline') + ax.plot(fl_spinup.dis_on_line, fl_spinup.bed_h, + label='Spinup Flowline') + ax.plot(fl_oggm.dis_on_line, fl_oggm.surface_h, + label='OGGM surface_h') + ax.plot(fl_spinup.dis_on_line, fl_spinup.surface_h, + label='Spinup surface_h') + # ax.plot(fl_consensus.dis_on_line, fl_consensus.surface_h, + # label='Consensus surface_h') + ax.plot(fl_agile_init.dis_on_line, fl_agile_init.surface_h, + label='Init surface_h') + ax.plot(fl_agile_end.dis_on_line, fl_agile_end.surface_h, + label='END surface_h') + ax.legend(); + + ax2.axhline(0, color='black') + ax2.plot(fl_spinup.surface_h - fl_agile_init.surface_h, + label='positive means retreating ' + '(spinup to rgi _date)') + ax2.plot(fl_agile_init.surface_h - fl_agile_end.surface_h, + label='rgi_date to end') + ax2.set_title('delta surface_h') + ax2.legend() + + ax3.axhline(0, color='black') + ax3.plot(fl_spinup.bed_h - fl_oggm_first_guess.bed_h, + label='positive means overdeepening') + ax3.set_title('delta bed_h') + ax3.legend() + + ax4.plot(fl_oggm_first_guess.is_trapezoid, label='trapez') + ax4.plot(fl_oggm_first_guess.is_rectangular, label='rect') + ax4.legend() + + ax5.plot(fl_oggm.thick > 0, label='OGGM thick > 0') + ax5.plot(fl_agile_init.thick > 0, + label='agile init thick > 0') + ax5.legend() + + ax6.plot(fl_diag_init.ice_velocity_myr[-1], + label=f'year {fl_diag_init.time[-1].values}') + ax6.plot(fl_diag_start.ice_velocity_myr[-1], + label=f'year {fl_diag_start.time[-1].values} ' + f'(rgi date)') + ax6.plot(fl_diag_end.ice_velocity_myr[-1], + label=f'year {fl_diag_end.time[-1].values}') + ax6.set_title('velocity') + ax6.legend() + + fig.suptitle(f'{gdir.name}, {glacier_state}', fontsize=16) + fig.tight_layout(pad=1.0) + plt.show() + + # test that glacier state mb is correctly used + inversion_settings = get_default_inversion_settings() + inversion_settings['minimize_options']['maxiter'] = 1 + + # add all possible observations to test everything + inversion_settings['observations']['fl_widths:m'] = {} + inversion_settings['obs_scaling_parameters']['uncertainty']['fl_widths:m'] = \ + {'absolute': 1.} + inversion_settings['observations']['fl_total_area:m2'] = {} + inversion_settings['obs_scaling_parameters']['uncertainty']['fl_total_area:m2'] = \ + {'absolute': 1.} + # extracted from experiment creation + inversion_settings['observations']['area:km2'] = {} + inversion_settings['obs_scaling_parameters']['uncertainty']['area:km2'] = \ + {'relative': 0.1} + + # add all possible control variables + inversion_settings['control_vars'] = ['area_bed_h'] + + for gdir in gdirs: + for glacier_state in glacier_states: + inversion_settings['experiment_description'] = ( + f"{glacier_state}_" + f"{inversion_settings['experiment_description']}") + data_logger = conduct_sandbox_inversion( + gdir, glacier_state=glacier_state, + inversion_settings=inversion_settings, + output_folder=test_dir, + init_model_fls='_oggm_first_guess', + gcm='BCC-CSM2-MR', + ssp='ssp370', + print_statistic=False, + return_data_logger=True) + + mb_models = initialise_mb_models(None, data_logger) + mb_model = mb_models['MB']['mb_model'].mbmod + mb_model_ref = get_glacier_state_mb_model(gdir, glacier_state) + mb_model_ref = mb_model_ref.flowline_mb_models[0] + + if glacier_state == 'equilibrium': + assert isinstance(mb_model, ConstantMassBalance) + assert mb_model.halfsize == mb_model_ref.halfsize, f"{gdir.rgi_id}" + assert mb_model.y0 == mb_model_ref.y0, f"{gdir.rgi_id}" + elif glacier_state in ['retreating', 'advancing']: + assert isinstance(mb_model, MonthlyTIModel) + assert mb_model.temp_bias == mb_model_ref.temp_bias, f"{gdir.rgi_id}" + assert mb_model.repeat == mb_model_ref.repeat, f"{gdir.rgi_id}" + assert mb_model.ys == mb_model_ref.ys, f"{gdir.rgi_id}" + assert mb_model.ye == mb_model_ref.ye, f"{gdir.rgi_id}" + else: + raise NotImplementedError(f'{glacier_state}') @pytest.mark.parametrize('control_vars', [['area_bed_h', 'lambdas', 'w0_m'], @@ -235,6 +317,7 @@ def test_run_idealized_experiment(self, test_dir, control_vars, init_mdl_fls): experiment_glacier = ['Aletsch', 'Artesonraju'] + glacier_states = ['equilibrium', 'retreating', 'advancing'] inversion_settings = get_default_inversion_settings() inversion_settings['minimize_options']['maxiter'] = 3 @@ -276,6 +359,7 @@ def test_run_idealized_experiment(self, test_dir, control_vars, gdirs = idealized_experiment( use_experiment_glaciers=experiment_glacier, + use_experiment_glacier_states=glacier_states, inversion_settings_all=[inversion_settings], inversion_settings_individual=inversion_settings_individual, init_model_fls=init_mdl_fls, @@ -286,93 +370,113 @@ def test_run_idealized_experiment(self, test_dir, control_vars, # test if individual inversion settings are correctly used for gdir in gdirs: - inv_setting = gdir.read_pickle(filename='inversion_input', - filesuffix='agile_inversion_results') - fg_hs = inv_setting['spinup_options']['height_shift']['mb_model']['fg_height_shift'] - if gdir.rgi_id == 'RGI60-11.01450': - assert fg_hs == -140 - elif gdir.rgi_id == 'RGI60-16.02444': - assert fg_hs == -52 - else: - raise NotImplementedError(f'{gdir.rgi_id}') - - # open final dataset - fp = os.path.join(test_dir, - 'Aletsch_agile_inversion_results.pkl') - with open(fp, 'rb') as handle: - ds = pickle.load(handle) + for glacier_state in glacier_states: + inv_setting = gdir.read_pickle( + filename='inversion_input', + filesuffix=f'_{glacier_state}_agile_inversion_results') + fg_hs = \ + inv_setting['spinup_options']['height_shift']['mb_model']['fg_height_shift'] + if gdir.rgi_id == 'RGI60-11.01450': + assert fg_hs == -140 + elif gdir.rgi_id == 'RGI60-16.02444': + assert fg_hs == -52 + else: + raise NotImplementedError(f'{gdir.rgi_id}') - # open default oggm statistics - fp = os.path.join(gdirs[0].dir, - 'default_oggm_statistics.pkl') - with open(fp, 'rb') as handle: - ds_default_stats = pickle.load(handle) - - # filesuffixes of two different oggm experiments - for oggm_run_suffix in ['_dynamic', '_static']: - # test for observation statistics - ds_key = 'observations_stats' - ds_key_oggm = ds_key + oggm_run_suffix - assert ds_key in ds.attrs.keys() - assert ds_key_oggm in ds_default_stats.keys() - for obs_key in ds.attrs[ds_key].keys(): - obs_key_name = obs_key.split(':')[0] - for year_key in ds.attrs[ds_key][obs_key].keys(): - if obs_key_name in ['fl_surface_h', 'fl_widths']: + for glacier_state in glacier_states: + # open final dataset + fp = os.path.join(test_dir, + f'Aletsch_{glacier_state}_agile_' + f'inversion_results.pkl') + with open(fp, 'rb') as handle: + ds = pickle.load(handle) + + # open default oggm statistics + fp = os.path.join(gdirs[0].dir, + f'default_oggm_statistics_{glacier_state}.pkl') + with open(fp, 'rb') as handle: + ds_default_stats = pickle.load(handle) + + # filesuffixes of two different oggm experiments + for oggm_run_suffix in ['_dynamic', '_static']: + # test for observation statistics + ds_key = 'observations_stats' + ds_key_oggm = ds_key + oggm_run_suffix + assert ds_key in ds.attrs.keys() + assert ds_key_oggm in ds_default_stats.keys() + for obs_key in ds.attrs[ds_key].keys(): + obs_key_name = obs_key.split(':')[0] + for year_key in ds.attrs[ds_key][obs_key].keys(): + if obs_key_name in ['fl_surface_h', 'fl_widths']: + test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] + elif obs_key_name in ['fl_total_area', 'area', 'dmdtda']: + test_metrics = ['diff', 'abs_diff'] + else: + raise NotImplementedError() + for metric in test_metrics: + assert metric in ds.attrs[ds_key][obs_key][year_key].keys() + assert metric in ds_default_stats[ds_key_oggm][obs_key][year_key].keys() + assert isinstance(ds.attrs[ds_key][obs_key][year_key][metric], + float) + assert isinstance(ds_default_stats[ds_key_oggm][obs_key][year_key][metric], + float) + + # test for control statistics + ds_key = 'controls_stats' + ds_key_oggm = ds_key + oggm_run_suffix + assert ds_key in ds.attrs.keys() + assert ds_key_oggm in ds_default_stats.keys() + for control_key in ds.attrs[ds_key].keys(): + if control_key in ['bed_h', 'area_bed_h', 'lambdas', 'w0_m']: test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] - elif obs_key_name in ['fl_total_area', 'area', 'dmdtda']: + test_default = True + elif control_key in ['height_shift_spinup']: test_metrics = ['diff', 'abs_diff'] + test_default = False else: raise NotImplementedError() for metric in test_metrics: - assert metric in ds.attrs[ds_key][obs_key][year_key].keys() - assert metric in ds_default_stats[ds_key_oggm][obs_key][year_key].keys() - assert isinstance(ds.attrs[ds_key][obs_key][year_key][metric], - float) - assert isinstance(ds_default_stats[ds_key_oggm][obs_key][year_key][metric], - float) - - # test for control statistics - ds_key = 'controls_stats' - ds_key_oggm = ds_key + oggm_run_suffix - assert ds_key in ds.attrs.keys() - assert ds_key_oggm in ds_default_stats.keys() - for control_key in ds.attrs[ds_key].keys(): - if control_key in ['bed_h', 'area_bed_h', 'lambdas', 'w0_m']: - test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] - test_default = True - elif control_key in ['height_shift_spinup']: - test_metrics = ['diff', 'abs_diff'] - test_default = False - else: - raise NotImplementedError() - for metric in test_metrics: - assert metric in ds.attrs[ds_key][control_key].keys() - assert isinstance(ds.attrs[ds_key][control_key][metric], - float) - if test_default: - assert metric in ds_default_stats[ds_key_oggm][control_key].keys() - assert isinstance(ds_default_stats[ds_key_oggm][control_key][metric], + assert metric in ds.attrs[ds_key][control_key].keys() + assert isinstance(ds.attrs[ds_key][control_key][metric], float) - - # test the past glacier state statistics - ds_key = 'past_state_stats' - ds_key_oggm = ds_key + oggm_run_suffix - assert ds_key in ds.attrs.keys() - assert ds_key_oggm in ds_default_stats.keys() - for var in ['thick', 'area_m2', 'volume_m3']: - test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'diff', 'bias'] - for metric in test_metrics: - if metric == 'diff': - if var != 'thick': + if test_default: + assert metric in ds_default_stats[ds_key_oggm][control_key].keys() + assert isinstance(ds_default_stats[ds_key_oggm][control_key][metric], + float) + + # test the past glacier state statistics + ds_key = 'past_state_stats' + ds_key_oggm = ds_key + oggm_run_suffix + assert ds_key in ds.attrs.keys() + assert ds_key_oggm in ds_default_stats.keys() + for var in ['thick', 'area_m2', 'volume_m3']: + test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'diff', 'bias'] + for metric in test_metrics: + if metric == 'diff': + if var != 'thick': + assert metric in ds.attrs[ds_key][var].keys() + assert metric in ds_default_stats[ds_key_oggm][var].keys() + assert isinstance(ds.attrs[ds_key][var][metric], + np.ndarray) + assert isinstance( + ds_default_stats[ds_key_oggm][var][metric], + np.ndarray) + else: assert metric in ds.attrs[ds_key][var].keys() assert metric in ds_default_stats[ds_key_oggm][var].keys() assert isinstance(ds.attrs[ds_key][var][metric], - np.ndarray) - assert isinstance( - ds_default_stats[ds_key_oggm][var][metric], - np.ndarray) - else: + float) + assert isinstance(ds_default_stats[ds_key_oggm][var][metric], + float) + + # test the past evolution statistics + ds_key = 'past_evol_stats' + ds_key_oggm = ds_key + oggm_run_suffix + assert ds_key in ds.attrs.keys() + assert ds_key_oggm in ds_default_stats.keys() + for var in ['volume_m3', 'area_m2']: + test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] + for metric in test_metrics: assert metric in ds.attrs[ds_key][var].keys() assert metric in ds_default_stats[ds_key_oggm][var].keys() assert isinstance(ds.attrs[ds_key][var][metric], @@ -380,38 +484,38 @@ def test_run_idealized_experiment(self, test_dir, control_vars, assert isinstance(ds_default_stats[ds_key_oggm][var][metric], float) - # test the past evolution statistics - ds_key = 'past_evol_stats' - ds_key_oggm = ds_key + oggm_run_suffix - assert ds_key in ds.attrs.keys() - assert ds_key_oggm in ds_default_stats.keys() - for var in ['volume_m3', 'area_m2']: - test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] - for metric in test_metrics: - assert metric in ds.attrs[ds_key][var].keys() - assert metric in ds_default_stats[ds_key_oggm][var].keys() - assert isinstance(ds.attrs[ds_key][var][metric], - float) - assert isinstance(ds_default_stats[ds_key_oggm][var][metric], - float) - - # test todays glacier state statistics - ds_key = 'today_state_stats' - ds_key_oggm = ds_key + oggm_run_suffix - assert ds_key in ds.attrs.keys() - assert ds_key_oggm in ds_default_stats.keys() - for var in ['thick', 'area_m2', 'volume_m3']: - test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'diff', 'bias'] - for metric in test_metrics: - if metric == 'diff': - if var != 'thick': + # test todays glacier state statistics + ds_key = 'today_state_stats' + ds_key_oggm = ds_key + oggm_run_suffix + assert ds_key in ds.attrs.keys() + assert ds_key_oggm in ds_default_stats.keys() + for var in ['thick', 'area_m2', 'volume_m3']: + test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'diff', 'bias'] + for metric in test_metrics: + if metric == 'diff': + if var != 'thick': + assert metric in ds.attrs[ds_key][var].keys() + assert metric in ds_default_stats[ds_key_oggm][var].keys() + assert isinstance(ds.attrs[ds_key][var][metric], + np.ndarray) + assert isinstance(ds_default_stats[ds_key_oggm][var][metric], + np.ndarray) + else: assert metric in ds.attrs[ds_key][var].keys() assert metric in ds_default_stats[ds_key_oggm][var].keys() assert isinstance(ds.attrs[ds_key][var][metric], - np.ndarray) + float) assert isinstance(ds_default_stats[ds_key_oggm][var][metric], - np.ndarray) - else: + float) + + # test the future evolution statistics + ds_key = 'future_evol_stats' + ds_key_oggm = ds_key + oggm_run_suffix + assert ds_key in ds.attrs.keys() + assert ds_key_oggm in ds_default_stats.keys() + for var in ['volume_m3', 'area_m2']: + test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] + for metric in test_metrics: assert metric in ds.attrs[ds_key][var].keys() assert metric in ds_default_stats[ds_key_oggm][var].keys() assert isinstance(ds.attrs[ds_key][var][metric], @@ -419,24 +523,10 @@ def test_run_idealized_experiment(self, test_dir, control_vars, assert isinstance(ds_default_stats[ds_key_oggm][var][metric], float) - # test the future evolution statistics - ds_key = 'future_evol_stats' - ds_key_oggm = ds_key + oggm_run_suffix - assert ds_key in ds.attrs.keys() - assert ds_key_oggm in ds_default_stats.keys() - for var in ['volume_m3', 'area_m2']: - test_metrics = ['rmsd', 'mean_ad', 'max_ad', 'bias'] - for metric in test_metrics: - assert metric in ds.attrs[ds_key][var].keys() - assert metric in ds_default_stats[ds_key_oggm][var].keys() - assert isinstance(ds.attrs[ds_key][var][metric], - float) - assert isinstance(ds_default_stats[ds_key_oggm][var][metric], - float) - def test_perfect_spinup_and_section_spinup(self, test_dir): experiment_glacier = ['Aletsch'] + glacier_states = ['retreating'] inversion_settings = get_default_inversion_settings() inversion_settings['minimize_options']['maxiter'] = 3 @@ -483,6 +573,7 @@ def test_perfect_spinup_and_section_spinup(self, test_dir): gdirs = idealized_experiment( use_experiment_glaciers=experiment_glacier, + use_experiment_glacier_states=glacier_states, inversion_settings_all=[inversion_settings1, inversion_settings2, inversion_settings3, @@ -494,10 +585,11 @@ def test_perfect_spinup_and_section_spinup(self, test_dir): 'cfl_number': 0.5}) fl_true_init = gdirs[0].read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] + filesuffix='_creation_spinup_' + f'{glacier_states[0]}')[0] # some tests for perfect sfc_h spinup fp = os.path.join(test_dir, - 'Aletsch_perfect_sfc_h_spinup.pkl') + f'Aletsch_{glacier_states[0]}_perfect_sfc_h_spinup.pkl') with open(fp, 'rb') as handle: ds_perfect_sfc_h = pickle.load(handle) @@ -510,7 +602,7 @@ def test_perfect_spinup_and_section_spinup(self, test_dir): # some tests for perfect thickness fp = os.path.join(test_dir, - 'Aletsch_perfect_thickness_spinup.pkl') + f'Aletsch_{glacier_states[0]}_perfect_thickness_spinup.pkl') with open(fp, 'rb') as handle: ds_perfect_thickness = pickle.load(handle) @@ -522,7 +614,7 @@ def test_perfect_spinup_and_section_spinup(self, test_dir): # some tests for perfect section fp = os.path.join(test_dir, - 'Aletsch_perfect_section_spinup.pkl') + f'Aletsch_{glacier_states[0]}_perfect_section_spinup.pkl') with open(fp, 'rb') as handle: ds_perfect_section = pickle.load(handle) @@ -538,12 +630,13 @@ def test_perfect_spinup_and_section_spinup(self, test_dir): # some tests for section spinup fp = os.path.join(test_dir, - 'Aletsch_section_spinup.pkl') + f'Aletsch_{glacier_states[0]}_section_spinup.pkl') with open(fp, 'rb') as handle: ds_section = pickle.load(handle) fl_fg = gdirs[0].read_pickle('model_flowlines', - filesuffix='_oggm_first_guess')[0] + filesuffix=f'_oggm_first_guess_' + f'{glacier_states[0]}')[0] assert np.allclose(ds_section.section_start[0], fl_fg.section) def test_StackedMassBalance(self, test_dir): @@ -561,10 +654,12 @@ def test_StackedMassBalance(self, test_dir): 'L1-L2_files/elev_bands/' glacier_names = ['Aletsch'] + glacier_states = ['retreating'] cfg.PARAMS['use_multiprocessing'] = False cfg.PARAMS['cfl_number'] = 0.5 gdirs = create_idealized_experiments(glacier_names, + glacier_states, prepro_border=prepro_border, from_prepro_level=from_prepro_level, base_url=base_url, ) @@ -624,6 +719,7 @@ def test_StackedMassBalance(self, test_dir): def test_perfect_bed_h(self, test_dir): experiment_glacier = ['Aletsch'] + glacier_states = ['retreating'] inversion_settings = get_default_inversion_settings() inversion_settings['minimize_options']['maxiter'] = 3 @@ -656,6 +752,7 @@ def test_perfect_bed_h(self, test_dir): gdirs = idealized_experiment( use_experiment_glaciers=experiment_glacier, + use_experiment_glacier_states=glacier_states, inversion_settings_all=[inversion_settings], working_dir=test_dir, output_folder=test_dir, @@ -663,10 +760,11 @@ def test_perfect_bed_h(self, test_dir): 'cfl_number': 0.5}) fl_true_init = gdirs[0].read_pickle('model_flowlines', - filesuffix='_creation_spinup')[0] + filesuffix='_creation_spinup_' + f'{glacier_states[0]}')[0] # some tests for perfect sfc_h spinup fp = os.path.join(test_dir, - 'Aletsch_perfect_bed_h.pkl') + f'Aletsch_{glacier_states[0]}_perfect_bed_h.pkl') with open(fp, 'rb') as handle: ds_perfect_bed_h = pickle.load(handle)