Initialise dynamic flowline with observed thickness

docs/whats-new.rst

@@ -6,10 +6,18 @@ Version history
 v1.6.2 (unreleased)
 -------------------
 
+Enhancements
+~~~~~~~~~~~~
+
+- There is now a possibility for initializing a elevation-band flowline using 
+  external thickness data and conduct a dynamic run with it (:pull:`1658`).
+
 Bug fixes
 ~~~~~~~~~
-- The binned variables in the elevation band flowlines did not used the
-  glacier mask when preserving the total values (:pull:`1661`).
+
+- The binned variables in the elevation band flowlines did not use the
+  glacier mask when preserving the total values. This is a bad 
+  bug that is now fixed (:pull:`1661`).
   By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 
 

oggm/core/flowline.py

@@ -2996,7 +2996,8 @@ def calving_glacier_downstream_line(line, n_points):
 
 
 @entity_task(log, writes=['model_flowlines'])
-def init_present_time_glacier(gdir, filesuffix=''):
+def init_present_time_glacier(gdir, filesuffix='',
+                              use_binned_thickness_data=False):
     """Merges data from preprocessing tasks. First task after inversion!
 
     This updates the `mode_flowlines` file and creates a stand-alone numerical
@@ -3007,9 +3008,14 @@ def init_present_time_glacier(gdir, filesuffix=''):
     gdir : :py:class:`oggm.GlacierDirectory`
         the glacier directory to process
     filesuffix : str
-            append a suffix to the model_flowlines filename (e.g. useful for
-            dynamic melt_f calibration including an inversion, so the original
-            model_flowlines are not changed).
+        append a suffix to the model_flowlines filename (e.g. useful for
+        dynamic melt_f calibration including an inversion, so the original
+        model_flowlines are not changed).
+    use_binned_thickness_data : bool or str
+        if you want to use thickness data, which was binned to the elevation
+        band flowlines with tasks.elevation_band_flowine and
+        tasks.fixed_dx_elevation_band_flowline, you can provide the name of the
+        data here to create a flowline for a dynamic model run
     """
 
     # Some vars
@@ -3026,20 +3032,53 @@ def init_present_time_glacier(gdir, filesuffix=''):
 
         # Get the data to make the model flowlines
         line = cl.line
-        section = inv['volume'] / (cl.dx * map_dx)
         surface_h = cl.surface_h
-        bed_h = surface_h - inv['thick']
         widths_m = cl.widths * map_dx
-
         assert np.all(widths_m > 0)
-        bed_shape = 4 * inv['thick'] / (cl.widths * map_dx) ** 2
 
-        lambdas = inv['thick'] * np.NaN
-        lambdas[inv['is_trapezoid']] = def_lambda
-        lambdas[inv['is_rectangular']] = 0.
+        if not use_binned_thickness_data:
+            # classical initialisation after the inversion
+            section = inv['volume'] / (cl.dx * map_dx)
+            bed_h = surface_h - inv['thick']
+
+            bed_shape = 4 * inv['thick'] / widths_m ** 2
+
+            lambdas = inv['thick'] * np.NaN
+            lambdas[inv['is_trapezoid']] = def_lambda
+            lambdas[inv['is_rectangular']] = 0.
 
-        # Where the flux and the thickness is zero we just assume trapezoid:
-        lambdas[bed_shape == 0] = def_lambda
+            # Where the flux and the thickness is zero we just assume trapezoid:
+            lambdas[bed_shape == 0] = def_lambda
+
+        else:
+            # here we use binned thickness data for the initialisation
+            elev_fl = pd.read_csv(
+                gdir.get_filepath('elevation_band_flowline',
+                                  filesuffix='_fixed_dx'), index_col=0)
+            assert np.allclose(widths_m, elev_fl['widths_m'])
+            elev_fl_thick = elev_fl[use_binned_thickness_data].values
+            section = elev_fl_thick * widths_m
+            lambdas = np.ones(len(section)) * def_lambda
+
+            # for trapezoidal the calculation of thickness results in quadratic
+            # equation, we only keep solution which results in a positive w0
+            # value (-> w/lambda >= h),
+            # it still could happen that we would need a negative w0 if the
+            # section is too large, in those cases we get a negative value
+            # inside sqrt (we ignore the RuntimeWarning) -> if this happens we
+            # use rectangular shape with original thickness
+            with np.errstate(invalid='ignore'):
+                thick = ((2 * widths_m -
+                          np.sqrt(4 * widths_m ** 2 -
+                                  4 * lambdas * 2 * section)) /
+                         (2 * lambdas))
+            nan_thick = np.isnan(thick)
+            thick[nan_thick] = elev_fl_thick[nan_thick]
+            lambdas[nan_thick] = 0
+
+            # finally the glacier bed and other stuff
+            bed_h = surface_h - thick
+            bed_shape = 4 * thick / widths_m ** 2
 
         if not gdir.is_tidewater and inv['is_last']:
             # for valley glaciers, simply add the downstream line, depending on

oggm/tests/test_models.py

@@ -17,7 +17,7 @@
 import xarray as xr
 from oggm import utils, workflow, tasks, cfg
 from oggm.core import climate, inversion, centerlines
-from oggm.shop import gcm_climate
+from oggm.shop import gcm_climate, bedtopo
 from oggm.cfg import SEC_IN_YEAR, SEC_IN_MONTH
 from oggm.utils import get_demo_file
 from oggm.exceptions import InvalidParamsError, InvalidWorkflowError
@@ -153,6 +153,67 @@ def test_init_present_time_glacier(self, hef_gdir, downstream_line_shape):
         with pytest.raises(InvalidParamsError):
             init_present_time_glacier(gdir)
 
+    def test_init_present_time_glacier_obs_thick(self, hef_elev_gdir,
+                                                 monkeypatch):
+
+        gdir = hef_elev_gdir
+
+        # need to change rgi_id, which is needed to be different in other tests
+        # when comparing to centerlines
+        gdir.rgi_id = 'RGI60-11.00897'
+
+        # add some thickness data
+        ft = utils.get_demo_file('RGI60-11.00897_thickness.tif')
+        monkeypatch.setattr(utils, 'file_downloader', lambda x: ft)
+        bedtopo.add_consensus_thickness(gdir)
+        vn = 'consensus_ice_thickness'
+        centerlines.elevation_band_flowline(gdir, bin_variables=[vn])
+        centerlines.fixed_dx_elevation_band_flowline(gdir,
+                                                     bin_variables=[vn])
+
+        tasks.init_present_time_glacier(gdir, filesuffix='_consensus',
+                                        use_binned_thickness_data=vn)
+        fl_consensus = gdir.read_pickle('model_flowlines',
+                                        filesuffix='_consensus')[0]
+
+        # check that resulting flowline has the same volume as observation
+        cdf = pd.read_hdf(utils.get_demo_file('rgi62_itmix_df.h5'))
+        ref_vol = cdf.loc[gdir.rgi_id].vol_itmix_m3
+        np.testing.assert_allclose(fl_consensus.volume_m3, ref_vol)
+
+        # should be trapezoid where ice
+        assert np.all(fl_consensus.is_trapezoid[fl_consensus.thick > 0])
+
+        # test that we can use fl in an dynamic model run without an error
+        mb = LinearMassBalance(3000.)
+        model_ref = FluxBasedModel(gdir.read_pickle('model_flowlines'),
+                                   mb_model=mb)
+        model_ref.run_until(100)
+        model_consensus = FluxBasedModel([fl_consensus], mb_model=mb)
+        model_consensus.run_until(100)
+        np.testing.assert_allclose(model_ref.volume_km3,
+                                   model_consensus.volume_km3,
+                                   atol=0.02)
+
+        # test that if w0<0 it is converted to rectangular
+        # set some thickness to very large values to force it
+        df_fixed_dx = pd.read_csv(gdir.get_filepath('elevation_band_flowline',
+                                                    filesuffix='_fixed_dx'))
+        new_thick = df_fixed_dx['consensus_ice_thickness']
+        new_thick[-10:] = new_thick[-10:] + 1000
+        df_fixed_dx['consensus_ice_thickness'] = new_thick
+        ref_vol_rect = np.sum(df_fixed_dx['area_m2'] * new_thick)
+        df_fixed_dx.to_csv(gdir.get_filepath('elevation_band_flowline',
+                                             filesuffix='_fixed_dx'))
+
+        tasks.init_present_time_glacier(gdir, filesuffix='_consensus_rect',
+                                        use_binned_thickness_data=vn)
+        fl_consensus_rect = gdir.read_pickle('model_flowlines',
+                                             filesuffix='_consensus_rect')[0]
+
+        np.testing.assert_allclose(fl_consensus_rect.volume_m3, ref_vol_rect)
+        assert np.sum(fl_consensus_rect.is_rectangular) == 10
+
     def test_present_time_glacier_massbalance(self, hef_gdir):
 
         gdir = hef_gdir