get ref mb profile with constant step size (#1366)

oggm/tests/test_models.py

@@ -148,16 +148,33 @@ def test_present_time_glacier_massbalance(self, hef_gdir):
         # Bias
         assert np.abs(utils.md(tot_mb, refmb)) < 50
 
+        mb_profile_constant_dh = gdir.get_ref_mb_profile(constant_dh=True)
+        step_dh = mb_profile_constant_dh.columns[1:] - mb_profile_constant_dh.columns[:-1]
+        assert np.all(step_dh == 50)
+        mb_profile_raw = gdir.get_ref_mb_profile()
+
         # Gradient
-        dfg = gdir.get_ref_mb_profile().mean()
+        dfg = mb_profile_raw.mean()
+        dfg_constant_dh = mb_profile_constant_dh.mean()
+
 
         # Take the altitudes below 3100 and fit a line
-        dfg = dfg[dfg.index < 3100]
         pok = np.where(hgts < 3100)
+        pok_constant_dh = np.where(dfg_constant_dh.index < 3100)
+
+        dfg = dfg[dfg.index < 3100]
+        dfg_constant_dh = dfg_constant_dh[dfg_constant_dh.index < 3100]
+
         from scipy.stats import linregress
         slope_obs, _, _, _, _ = linregress(dfg.index, dfg.values)
+        slope_obs_constant_dh, _, _, _, _ = linregress(dfg_constant_dh.index,
+                                                       dfg_constant_dh.values)
         slope_our, _, _, _, _ = linregress(hgts[pok], grads[pok])
+
         np.testing.assert_allclose(slope_obs, slope_our, rtol=0.15)
+        # the observed MB gradient and the interpolated observed one (with
+        # constant dh) should be very similar!
+        np.testing.assert_allclose(slope_obs, slope_obs_constant_dh, rtol=0.01)
 
 
 @pytest.fixture(scope='class')

oggm/utils/_workflow.py

@@ -2482,6 +2482,7 @@ def __init__(self, rgi_entity, base_dir=None, reset=False,
         # Optimization
         self._mbdf = None
         self._mbprofdf = None
+        self._mbprofdf_cte_dh = None
 
     def __repr__(self):
 
@@ -3216,7 +3217,7 @@ def get_ref_mb_data(self, y0=None, y1=None, input_filesuffix=''):
             out = self._mbdf
         return out.dropna(subset=['ANNUAL_BALANCE'])
 
-    def get_ref_mb_profile(self, input_filesuffix=''):
+    def get_ref_mb_profile(self, input_filesuffix='', constant_dh=False):
         """Get the reference mb profile data from WGMS (if available!).
 
         Returns None if this glacier has no profile and an Error if it isn't
@@ -3227,10 +3228,14 @@ def get_ref_mb_profile(self, input_filesuffix=''):
         input_filesuffix : str
             input_filesuffix of the climate_historical that should be used. The
             default is to take the climate_historical without input_filesuffix
-
+        constant_dh : boolean
+            If set to True, it outputs the MB profiles with a constant step size
+            of dh=50m by using interpolation. This can be useful for comparisons
+            between years. Default is False which gives the raw
+            elevation-dependent point MB
         """
 
-        if self._mbprofdf is None:
+        if self._mbprofdf is None and not constant_dh:
             flink, mbdatadir = get_wgms_files()
             c = 'RGI{}0_ID'.format(self.rgi_version[0])
             wid = flink.loc[flink[c] == self.rgi_id]
@@ -3247,16 +3252,40 @@ def get_ref_mb_profile(self, input_filesuffix=''):
             # list of years
             self._mbprofdf = pd.read_csv(reff, index_col=0)
 
+        if self._mbprofdf_cte_dh is None and constant_dh:
+            flink, mbdatadir = get_wgms_files()
+            c = 'RGI{}0_ID'.format(self.rgi_version[0])
+            wid = flink.loc[flink[c] == self.rgi_id]
+            if len(wid) == 0:
+                raise RuntimeError('Not a reference glacier!')
+            wid = wid.WGMS_ID.values[0]
+
+            # file
+            mbdatadir = os.path.join(os.path.dirname(mbdatadir), 'mb_profiles_constant_dh')
+            reff = os.path.join(mbdatadir,
+                                'profile_constant_dh_WGMS-{:05d}.csv'.format(wid))
+            if not os.path.exists(reff):
+                return None
+            # list of years
+            self._mbprofdf_cte_dh = pd.read_csv(reff, index_col=0)
+
         ci = self.get_climate_info(input_filesuffix=input_filesuffix)
         if 'baseline_hydro_yr_0' not in ci:
             raise RuntimeError('Please process some climate data before call')
         y0 = ci['baseline_hydro_yr_0']
         y1 = ci['baseline_hydro_yr_1']
-        if len(self._mbprofdf) > 1:
-            out = self._mbprofdf.loc[y0:y1]
+        if not constant_dh:
+            if len(self._mbprofdf) > 1:
+                out = self._mbprofdf.loc[y0:y1]
+            else:
+                # Some files are just empty
+                out = self._mbprofdf
         else:
-            # Some files are just empty
-            out = self._mbprofdf
+            if len(self._mbprofdf_cte_dh) > 1:
+                out = self._mbprofdf_cte_dh.loc[y0:y1]
+            else:
+                # Some files are just empty
+                out = self._mbprofdf_cte_dh
         out.columns = [float(c) for c in out.columns]
         return out.dropna(axis=1, how='all').dropna(axis=0, how='all')