From 59d300875e2304eb8f78403c8586e2da6ed824e3 Mon Sep 17 00:00:00 2001 From: Patrick Date: Fri, 21 Jul 2023 08:35:17 +0200 Subject: [PATCH] added second option for a worse first guess --- agile1d/core/inversion.py | 10 ++-- agile1d/sandbox/calculate_statistics.py | 2 +- .../create_glaciers_with_measurements.py | 51 +++++++++++++++++-- .../sandbox/define_idealized_experiment.py | 32 ++++++------ agile1d/sandbox/graphics.py | 2 +- agile1d/sandbox/run_idealized_experiment.py | 11 ++++ agile1d/tests/test_sandbox.py | 24 ++++++--- 7 files changed, 100 insertions(+), 32 deletions(-) diff --git a/agile1d/core/inversion.py b/agile1d/core/inversion.py index 1e481e3..8d606db 100644 --- a/agile1d/core/inversion.py +++ b/agile1d/core/inversion.py @@ -490,11 +490,11 @@ def save_past_evolution_to_disk(gdir, data_logger): @entity_task(log, writes=['model_flowlines']) def agile_inversion(gdir, inversion_input_filesuffix='_agile', - init_model_filesuffix=None, init_model_fls='_trapezoidal', - climate_filename='climate_historical', - climate_filesuffix='', output_filesuffix='_agile_output', - output_filepath=None, save_dataset=True, - give_data_logger_back=False, save_past_evolution=True): + init_model_filesuffix=None, init_model_fls='_trapezoidal', + climate_filename='climate_historical', + climate_filesuffix='', output_filesuffix='_agile_output', + output_filepath=None, save_dataset=True, + give_data_logger_back=False, save_past_evolution=True): """TODO """ log.debug('initialise Datalogger') diff --git a/agile1d/sandbox/calculate_statistics.py b/agile1d/sandbox/calculate_statistics.py index 6aa3682..bfaf9eb 100644 --- a/agile1d/sandbox/calculate_statistics.py +++ b/agile1d/sandbox/calculate_statistics.py @@ -373,7 +373,7 @@ def calculate_default_oggm_statistics(gdir): controls_mdl = {} controls_true = {} fls_mdl = gdir.read_pickle('model_flowlines', - filesuffix='_agile_first_guess')[0] + filesuffix='_oggm_first_guess')[0] fls_true = gdir.read_pickle('model_flowlines', filesuffix='_agile_true_init')[0] diff --git a/agile1d/sandbox/create_glaciers_with_measurements.py b/agile1d/sandbox/create_glaciers_with_measurements.py index c672798..ac5f82e 100644 --- a/agile1d/sandbox/create_glaciers_with_measurements.py +++ b/agile1d/sandbox/create_glaciers_with_measurements.py @@ -12,6 +12,7 @@ from oggm.shop import bedtopo, gcm_climate from agile1d.core.massbalance import StackedMassBalance +from agile1d.core.inversion import get_adaptive_upper_ice_thickness_limit from agile1d.sandbox.calculate_statistics import calculate_default_oggm_statistics from agile1d.sandbox.glaciers_for_idealized_experiments import experiment_glaciers @@ -131,6 +132,9 @@ def create_idealized_experiments(all_glaciers, # 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, @@ -538,7 +542,48 @@ 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='_agile_first_guess') + filesuffix='_oggm_first_guess') + + +@entity_task(log, writes=['model_flowlines']) +def glabtop_inversion_for_first_guess(gdir): + """TODO + """ + fls_experiment = gdir.read_pickle('model_flowlines', + filesuffix='_agile_true_init')[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) + fg_bed_h = fls_experiment.surface_h[fls_experiment.thick > 0] - fg_thick + + bed_h_new = copy.deepcopy(fls_experiment.bed_h) + bed_h_new[ice_mask] = fg_bed_h + + section_new = np.zeros(fls_experiment.nx) + widths = fls_experiment.widths_m[ice_mask] + lam = fls_experiment._lambdas[ice_mask] + section_new[ice_mask] = (2*widths - lam * fg_thick) / 2 * fg_thick + + fl_fg = MixedBedFlowline(line=fls_experiment.line, + dx=fls_experiment.dx, + map_dx=fls_experiment.map_dx, + surface_h=fls_experiment.surface_h, + bed_h=bed_h_new, + section=section_new, + bed_shape=fls_experiment.bed_shape, + is_trapezoid=fls_experiment.is_trapezoid, + lambdas=fls_experiment._lambdas, + widths_m=fls_experiment.widths_m, + rgi_id=fls_experiment.rgi_id, + 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.is_trapezoid) + + gdir.write_pickle([fl_fg], 'model_flowlines', + filesuffix='_glabtop_first_guess') def get_experiment_mb_model(gdir): @@ -556,7 +601,7 @@ def oggm_default_initialisation(gdir, ys, ye, gcm='BCC-CSM2-MR', ssp='ssp370'): # do dynamic initialisation spinup_model = tasks.run_dynamic_spinup( gdir, spinup_start_yr=ys, ye=ye, evolution_model=SemiImplicitModel, - model_flowline_filesuffix='_agile_first_guess', + model_flowline_filesuffix='_oggm_first_guess', precision_absolute=0.1, output_filesuffix='_oggm_dynamic_past', store_model_geometry=True, store_fl_diagnostics=True, store_model_evolution=True, @@ -564,7 +609,7 @@ def oggm_default_initialisation(gdir, ys, ye, gcm='BCC-CSM2-MR', ssp='ssp370'): # do static initialisation fls_init = gdir.read_pickle('model_flowlines', - filesuffix='_agile_first_guess') + filesuffix='_oggm_first_guess') static_model = tasks.run_from_climate_data( gdir, ye=ye, fixed_geometry_spinup_yr=ys, diff --git a/agile1d/sandbox/define_idealized_experiment.py b/agile1d/sandbox/define_idealized_experiment.py index cb99cd6..eadf522 100644 --- a/agile1d/sandbox/define_idealized_experiment.py +++ b/agile1d/sandbox/define_idealized_experiment.py @@ -21,6 +21,7 @@ def idealized_experiment(use_experiment_glaciers, working_dir, output_folder, inversion_settings_individual=None, params_file=None, override_params=None, + init_model_fls='_oggm_first_guess', logging_level='WORKFLOW', gcm='BCC-CSM2-MR', ssp='ssp370', @@ -81,6 +82,7 @@ def idealized_experiment(use_experiment_glaciers, )) workflow.execute_entity_task(conduct_sandbox_inversion, all_experiments, + init_model_fls=init_model_fls, gcm=gcm, ssp=ssp, print_statistic=print_statistic) @@ -116,7 +118,7 @@ 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, output_folder=None, - init_model_fls='_agile_first_guess', + init_model_fls='_oggm_first_guess', gcm='BCC-CSM2-MR', ssp='ssp370', print_statistic=True): @@ -143,22 +145,22 @@ def conduct_sandbox_inversion(gdir, inversion_settings=None, raise RuntimeError(f'No observation for {msr} available!') prepare_for_agile_inversion(gdir, - inversion_settings=inversion_settings, - filesuffix= - inversion_settings['experiment_description']) + inversion_settings=inversion_settings, + filesuffix= + inversion_settings['experiment_description']) data_logger = agile_inversion(gdir, - inversion_input_filesuffix= - inversion_settings['experiment_description'], - 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_filepath=output_folder, - save_dataset=True, - give_data_logger_back=True) + inversion_input_filesuffix= + inversion_settings['experiment_description'], + 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_filepath=output_folder, + save_dataset=True, + give_data_logger_back=True) add_future_projection_run(gdir, data_logger=data_logger, gcm=gcm, ssp=ssp) diff --git a/agile1d/sandbox/graphics.py b/agile1d/sandbox/graphics.py index eb8bf4d..732682c 100644 --- a/agile1d/sandbox/graphics.py +++ b/agile1d/sandbox/graphics.py @@ -503,7 +503,7 @@ def plot_glacier_evolution(gdir, ax_v, ax_a, ax_l, text_position=1.1, def plot_dhdt_first_guess_db(gdir, 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='_agile_first_guess')[0] + fls_first_guess = gdir.read_pickle('model_flowlines', filesuffix='_oggm_first_guess')[0] # define extend to plot according to longest flowline diff --git a/agile1d/sandbox/run_idealized_experiment.py b/agile1d/sandbox/run_idealized_experiment.py index 70fe703..44584bc 100644 --- a/agile1d/sandbox/run_idealized_experiment.py +++ b/agile1d/sandbox/run_idealized_experiment.py @@ -35,6 +35,8 @@ def parse_args(args): default=False, help='If the idealized statistics should be printed ' 'out after each run.') + parser.add_argument('--first_guess', type=str, default='oggm', + help='The first guess method to use (oggm, glabtop).') args = parser.parse_args(args) @@ -49,6 +51,14 @@ def parse_args(args): output_folder = os.path.abspath(output_folder) working_dir = os.path.abspath(working_dir) + # define first guess flowline to use + if args.first_guess == 'oggm': + init_model_fls = '_oggm_first_guess' + elif args.first_guess == 'glabtop': + init_model_fls = '_glabtop_first_guess' + else: + raise NotImplementedError(f'{args.first_guess}') + # check if information for experiments is correctly given in separate file spec = importlib.util.spec_from_file_location("experiments", args.experiment_file) @@ -74,6 +84,7 @@ def parse_args(args): logging_level=args.logging_level, inversion_settings_all=foo.inversion_settings_all, inversion_settings_individual=foo.inversion_settings_individual, + init_model_fls=init_model_fls, use_experiment_glaciers=use_experiment_glaciers, print_statistic=args.print_statistic ) diff --git a/agile1d/tests/test_sandbox.py b/agile1d/tests/test_sandbox.py index 3bc1b16..463ecbb 100644 --- a/agile1d/tests/test_sandbox.py +++ b/agile1d/tests/test_sandbox.py @@ -135,6 +135,11 @@ def all_stats_finite(ds, use_year, is_static=False): 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') + if do_plot: for gdir in gdirs: fl_oggm = gdir.read_pickle('model_flowlines')[0] @@ -148,9 +153,9 @@ def all_stats_finite(ds, use_year, is_static=False): fl_agile_end = \ gdir.read_pickle('model_flowlines', filesuffix='_agile_true_end')[0] - fl_agile_first_guess = \ + fl_oggm_first_guess = \ gdir.read_pickle('model_flowlines', - filesuffix='_agile_first_guess')[0] + filesuffix='_oggm_first_guess')[0] def get_fl_diagnostics(filesuffix): f = gdir.get_filepath('fl_diagnostics', @@ -191,13 +196,13 @@ def get_fl_diagnostics(filesuffix): ax2.legend() ax3.axhline(0, color='black') - ax3.plot(fl_spinup.bed_h - fl_agile_first_guess.bed_h, + 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_agile_first_guess.is_trapezoid, label='trapez') - ax4.plot(fl_agile_first_guess.is_rectangular, label='rect') + 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') @@ -223,7 +228,11 @@ def get_fl_diagnostics(filesuffix): [['area_bed_h', 'lambdas', 'w0_m'], ['bed_h', 'lambdas', 'w0_m']], ids=['area_bed_h', 'bed_h']) - def test_run_idealized_experiment(self, test_dir, control_vars): + @pytest.mark.parametrize('init_mdl_fls', + ['_oggm_first_guess', '_glabtop_first_guess'], + ids=['oggm_fg', 'glabtop_fg']) + def test_run_idealized_experiment(self, test_dir, control_vars, + init_mdl_fls): experiment_glacier = ['Aletsch', 'Artesonraju'] @@ -269,6 +278,7 @@ def test_run_idealized_experiment(self, test_dir, control_vars): use_experiment_glaciers=experiment_glacier, inversion_settings_all=[inversion_settings], inversion_settings_individual=inversion_settings_individual, + init_model_fls=init_mdl_fls, working_dir=test_dir, output_folder=test_dir, override_params={'border': 160, @@ -533,7 +543,7 @@ def test_perfect_spinup_and_section_spinup(self, test_dir): ds_section = pickle.load(handle) fl_fg = gdirs[0].read_pickle('model_flowlines', - filesuffix='_agile_first_guess')[0] + filesuffix='_oggm_first_guess')[0] assert np.allclose(ds_section.section_start[0], fl_fg.section) def test_StackedMassBalance(self, test_dir):