Allow for monthly temp and precip biases / facs (#1247)

* WIP

* Allow for monthly temp and precip biases / facs

* Whats new and better handling of scalars

Author: fmaussion

docs/whats-new.rst

@@ -4,7 +4,7 @@
 Version history
 ===============
 
-v1.4.1 (unreleased)
+v1.5.0 (unreleased)
 -------------------
 
 Breaking changes
@@ -44,6 +44,9 @@ Enhancements
 - Added new diagnostic variables (`terminus_thick_0` ...) to better track
   the thickness at the terminus (:pull:`1230`).
   By `Fabien Maussion <https://github.com/fmaussion>`_
+- Temperature and precipitation corrections can now also be applied on a
+  monthly basis for the PastMassBalanceModel (:pull:`1247`).
+  By `Fabien Maussion <https://github.com/fmaussion>`_
 
 
 Bug fixes
@@ -52,8 +55,8 @@ Bug fixes
 - Fixed bug in hydro date / calendar date conversions with month=1 (i.e.
   no conversion) (:pull:`1220`).
   By `Lilian Schuster <https://github.com/lilianschuster>`_
-- Fixed bug in ``graphics.plot_distributed_thickness`` which encounter Error under
-  the elevation bands flowline (:pull: `1241`).
+- Fixed bug in ``graphics.plot_distributed_thickness`` which led to an error with
+  elevation bands flowlines (:pull:`1241`).
   By `Li Fei <https://github.com/Keeptg>`_
 
 v1.4.0 (17.02.2021)

oggm/core/massbalance.py

@@ -428,9 +428,18 @@ def prcp_fac(self):
     @prcp_fac.setter
     def prcp_fac(self, new_prcp_fac):
         # just to check that no invalid prcp_factors are used
-        if new_prcp_fac <= 0:
+        if np.any(np.asarray(new_prcp_fac) <= 0):
             raise InvalidParamsError('prcp_fac has to be above zero!')
+
+        if len(np.atleast_1d(new_prcp_fac)) == 12:
+            # OK so that's monthly stuff
+            # We dirtily assume that user just used calendar month
+            sm = cfg.PARAMS['hydro_month_' + self.hemisphere]
+            new_prcp_fac = np.roll(new_prcp_fac, 13 - sm)
+            new_prcp_fac = np.tile(new_prcp_fac, len(self.prcp) // 12)
+
         self.prcp *= new_prcp_fac / self._prcp_fac
+
         # update old prcp_fac in order that it can be updated again ...
         self._prcp_fac = new_prcp_fac
 
@@ -441,7 +450,16 @@ def temp_bias(self):
 
     @temp_bias.setter
     def temp_bias(self, new_temp_bias):
+
+        if len(np.atleast_1d(new_temp_bias)) == 12:
+            # OK so that's monthly stuff
+            # We dirtily assume that user just used calendar month
+            sm = cfg.PARAMS['hydro_month_' + self.hemisphere]
+            new_temp_bias = np.roll(new_temp_bias, 13 - sm)
+            new_temp_bias = np.tile(new_temp_bias, len(self.temp) // 12)
+
         self.temp += new_temp_bias - self._temp_bias
+
         # update old temp_bias in order that it can be updated again ...
         self._temp_bias = new_temp_bias
 

oggm/tests/test_models.py

@@ -393,6 +393,12 @@ def test_prcp_fac_temp_bias_update(self, hef_gdir):
         assert mb_mod._prcp_fac == prcp_fac_old
         assert mb_mod.temp_bias == temp_bias_old
 
+        # Now monthly stuff
+        mb_mod.temp_bias = [0] * 12
+        np.testing.assert_allclose(mb_mod.temp_bias, temp_bias_old)
+        mb_mod.prcp_fac = [prcp_fac_old] * 12
+        np.testing.assert_allclose(mb_mod.prcp_fac, prcp_fac_old)
+
         # Deprecated attrs
         with pytest.raises(AttributeError):
             mb_mod.prcp_bias = 2
@@ -560,6 +566,7 @@ def test_constant_mb_model(self, hef_gdir):
         months = np.arange(12)
         monthly_1 = months * 0.
         monthly_2 = months * 0.
+        monthly_3 = months * 0.
         for m in months:
             yr = utils.date_to_floatyear(0, m + 1)
             cmb_mod.temp_bias = 0
@@ -568,14 +575,23 @@ def test_constant_mb_model(self, hef_gdir):
             cmb_mod.temp_bias = 1
             tmp = cmb_mod.get_monthly_mb(h, yr) * SEC_IN_MONTH * rho
             monthly_2[m] = np.average(tmp, weights=w)
+            cmb_mod.temp_bias = [0] * 6 + [1] + [0] * 5
+            cmb_mod.prcp_fac = [10] * 3 + [2.5] * 9
+            tmp = cmb_mod.get_monthly_mb(h, yr) * SEC_IN_MONTH * rho
+            monthly_3[m] = np.average(tmp, weights=w)
+            cmb_mod.prcp_fac = 2.5
 
         # check that the winter months are close but summer months no
         np.testing.assert_allclose(monthly_1[1: 5], monthly_2[1: 5], atol=1)
+        np.testing.assert_allclose(monthly_1[1: 2], monthly_3[1: 2], atol=1)
         assert np.mean(monthly_1[5:]) > (np.mean(monthly_2[5:]) + 100)
+        assert np.mean(monthly_3[3:5]) > (np.mean(monthly_1[3:5]) + 100)
+        assert monthly_3[9] == monthly_2[9]
 
         if do_plot:  # pragma: no cover
             plt.plot(monthly_1, '-', label='Normal')
             plt.plot(monthly_2, '-', label='Temp bias')
+            plt.plot(monthly_3, '-', label='Temp bias monthly')
             plt.legend()
             plt.show()