added option to use gridded 2D topo data for geodetic_mb_calibration (#…

…1709)

* added usage of 2D data to mb_calibration_from_geodetic_mb

* added a test case for the new 2D mass calibration option

* PEP8 corrections

* added fabi's adaptions

* Updated whats-new.rst

---------

Co-authored-by: afisc <[email protected]>

docs/whats-new.rst

@@ -33,6 +33,11 @@ Enhancements
   Based on :pull:`1403` and :pull:`1532` (unmerged).
   By `Jan Malles <https://github.com/jmalles>`_ and
   `Fabien Maussion <https://github.com/fmaussion>`_
+- Added the option to use the 2D mask of a glacier for the calibration 
+  of the mass balance. One can now calibrate the mass balance from 2D data, which 
+  adds consistency when also the run itself is fully distributed. Changes are 
+  found in :pull:`1709`.
+  By `Alex Fischer <https://github.com/afisc>`_
 
 Bug fixes
 ~~~~~~~~~

oggm/core/massbalance.py

@@ -1433,6 +1433,7 @@ def mb_calibration_from_geodetic_mb(gdir, *,
                                     overwrite_gdir=False,
                                     use_regional_avg=False,
                                     override_missing=None,
+                                    use_2d_mb=False,
                                     informed_threestep=False,
                                     calibrate_param1='melt_f',
                                     calibrate_param2=None,
@@ -1479,6 +1480,9 @@ def mb_calibration_from_geodetic_mb(gdir, *,
         if the reference geodetic data is not available, use this value instead
         (mostly for testing with exotic datasets, but could be used to open
         the door to using other datasets).
+    use_2d_mb : bool
+        Set to True if the mass balance calibration has to be done of the 2D mask
+        of the glacier (for fully distributed runs only).
     informed_threestep : bool
         the magic method Fabi found out one day before release.
         Overrides the calibrate_param order below.
@@ -1574,6 +1578,7 @@ def mb_calibration_from_geodetic_mb(gdir, *,
                                              ref_period=ref_period,
                                              write_to_gdir=write_to_gdir,
                                              overwrite_gdir=overwrite_gdir,
+                                             use_2d_mb=use_2d_mb,
                                              calibrate_param1='prcp_fac',
                                              calibrate_param2='melt_f',
                                              calibrate_param3='temp_bias',
@@ -1592,6 +1597,7 @@ def mb_calibration_from_geodetic_mb(gdir, *,
                                              ref_period=ref_period,
                                              write_to_gdir=write_to_gdir,
                                              overwrite_gdir=overwrite_gdir,
+                                             use_2d_mb=use_2d_mb,
                                              calibrate_param1=calibrate_param1,
                                              calibrate_param2=calibrate_param2,
                                              calibrate_param3=calibrate_param3,
@@ -1609,6 +1615,7 @@ def mb_calibration_from_scalar_mb(gdir, *,
                                   ref_mb_years=None,
                                   write_to_gdir=True,
                                   overwrite_gdir=False,
+                                  use_2d_mb=False,
                                   calibrate_param1='melt_f',
                                   calibrate_param2=None,
                                   calibrate_param3=None,
@@ -1673,6 +1680,9 @@ def mb_calibration_from_scalar_mb(gdir, *,
         if a `mb_calib.json` exists, this task won't overwrite it per default.
         Set this to True to enforce overwriting (i.e. with consequences for the
         future workflow).
+    use_2d_mb : bool
+        Set to True if the mass balance calibration has to be done of the 2D mask
+        of the glacier (for fully distributed runs only).
     mb_model_class : MassBalanceModel class
         the MassBalanceModel to use for the calibration. Needs to use the
         same parameters as MonthlyTIModel (the default): melt_f,
@@ -1741,7 +1751,17 @@ def mb_calibration_from_scalar_mb(gdir, *,
         raise InvalidParamsError('Cannot set `ref_mb_years` and `ref_period` '
                                  'at the same time.')
 
-    fls = gdir.read_pickle('inversion_flowlines')
+    if not use_2d_mb:
+        fls = gdir.read_pickle('inversion_flowlines')
+    else:
+        # if the 2D data is used, the flowline is not needed.
+        fls = None
+        # get the 2D data
+        fp = gdir.get_filepath('gridded_data')
+        with xr.open_dataset(fp) as ds:
+            # 'topo' instead of 'topo_smoothed'?
+            heights = ds.topo_smoothed.data[ds.glacier_mask.data == 1]
+            widths = np.ones(len(heights))
 
     # Let's go
     # Climate period
@@ -1814,7 +1834,10 @@ def mb_calibration_from_scalar_mb(gdir, *,
     def to_minimize(x, model_attr):
         # Set the new attr value
         setattr(mb_mod, model_attr, x)
-        out = mb_mod.get_specific_mb(fls=fls, year=years).mean()
+        if use_2d_mb:
+            out = mb_mod.get_specific_mb(heights=heights, widths=widths, year=years).mean()
+        else:
+            out = mb_mod.get_specific_mb(fls=fls, year=years).mean()
         return np.mean(out - ref_mb)
 
     try:

oggm/tests/test_prepro.py

@@ -1699,6 +1699,21 @@ def test_mb_calibration_from_scalar_mb(self):
         assert pdf['melt_f'] < cfg.PARAMS['melt_f']
         assert pdf['prcp_fac'] == cfg.PARAMS['prcp_fac_max']
 
+        # Test the use of gridded data(2D) instead of flowline data(1D) for the calibration
+        mb_calibration_from_scalar_mb(gdir,
+                                      ref_mb=ref_mb,
+                                      ref_period=ref_period,
+                                      use_2d_mb=False)
+        mb_calib_1d = gdir.read_json('mb_calib')
+
+        mb_calibration_from_scalar_mb(gdir,
+                                      ref_mb=ref_mb,
+                                      ref_period=ref_period,
+                                      use_2d_mb=True)
+        mb_calib_2d = gdir.read_json('mb_calib')
+        # the calibration results for the melt factor should be close to each other (+/- 5% are tolerated)
+        np.testing.assert_allclose(mb_calib_2d['melt_f'], mb_calib_1d['melt_f'], rtol=0.05)
+
     @pytest.mark.slow
     def test_mb_calibration_from_scalar_mb_multiple_fl(self):