added second option for a worse first guess

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')

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]
 

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,15 +601,15 @@ 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,
         ignore_errors=False, add_fixed_geometry_spinup=True)
 
     # 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,

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)
 

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

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
                 )

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):