diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 4d1ac8603..3ffd4202c 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1445,7 +1445,7 @@ def mb_calibration_from_scalar_mb(gdir, calibrate_param1='melt_f', calibrate_param2=None, calibrate_param3=None, - melt_f_default=None, + melt_f=None, melt_f_min=None, melt_f_max=None, prcp_scaling_factor=None, @@ -1458,9 +1458,27 @@ def mb_calibration_from_scalar_mb(gdir, ): """Determine the mass balance parameters from a scalar mass-balance value. + If no reference value is given, we will use data from Hugonnet et al., 2021. + The data table can be obtained with utils.get_geodetic_mb_dataframe(). + It is equivalent to the original data, but has some outlier values + filtered. See `this notebook <>`_ for more details. + This calibrates the mass balance parameters using a reference average MB data over a given period (e.g. average in-situ SMB or geodetic MB). - This flexible calibration allows to calibrate three parameters. + This flexible calibration allows to calibrate three parameters one after + another. The first parameter is varied between two chosen values (a range) + until the ref MB value is matched. If this fails, the second parameter + can be changed, etc. + + This is similar to the "three-step calibration" introduced by + Huss & Hock 2015, but with some differences: + - this method is flexible, i.e. you can choose any order of calibration + - we use different defaults for the temperature bias and the + precipitation factor (see documentation) + + As a result, the current default in OGGM is to chose certain defaults + for prcp_fac and temp_bias, and calibrate melt_f. If this fails, + calibrate temp_bias. Note that this does not compute the apparent mass balance at the same time - users need to run `apparent_mb_from_any_mb after` @@ -1492,22 +1510,54 @@ def mb_calibration_from_scalar_mb(gdir, Set this to True to enforce overwriting (i.e. with consequences for the future workflow). override_missing : scalar - if the reference geodetic data is not available - calibrate_param1='melt_f', - calibrate_param2=None, - calibrate_param3=None, - melt_f_default=None, - melt_f_min=None, - melt_f_max=None, - prcp_scaling_factor=None, - prcp_scaling_factor_min=None, - prcp_scaling_factor_max=None, - temp_bias=0, - temp_bias_min=None, - temp_bias_max=None, + 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). mb_model_class : MassBalanceModel class - the MassBalanceModel to use for the calibration - default is MonthlyTIModel + the MassBalanceModel to use for the calibration. Needs to use the + same parameters as MonthlyTIModel (the default): melt_f, + temp_bias, prcp_fac. + calibrate_param1='melt_f' : str + in the three-step calibration, the name of the first parameter + to calibrate (one of 'melt_f', 'temp_bias', 'prcp_fac'). + calibrate_param2 : str + in the three-step calibration, the name of the second parameter + to calibrate (one of 'melt_f', 'temp_bias', 'prcp_fac'). If not + set and the algorithm cannot match observations, it will raise an + error. + calibrate_param3 : str + in the three-step calibration, the name of the third parameter + to calibrate (one of 'melt_f', 'temp_bias', 'prcp_fac'). If not + set and the algorithm cannot match observations, it will raise an + error. + melt_f: float + the default value to use as melt factor (or the starting value when + optimizing MB). Defaults to cfg.PARAMS['melt_f']. + melt_f_min: float + the minimum accepted value for the melt factor during optimisation. + Defaults to cfg.PARAMS['melt_f_min']. + melt_f_max: float + the maximum accepted value for the melt factor during optimisation. + Defaults to cfg.PARAMS['melt_f_max']. + prcp_scaling_factor: float + the default value to use as precipitation scaling factor + (or the starting value when optimizing MB). Defaults to the method + chosen in `params.cfg` (winter prcp or global factor). + prcp_scaling_factor_min: float + the minimum accepted value for the precipitation scaling factor during + optimisation. Defaults to cfg.PARAMS['prcp_scaling_factor_min']. + prcp_scaling_factor_max: float + the maximum accepted value for the precipitation scaling factor during + optimisation. Defaults to cfg.PARAMS['prcp_scaling_factor_max']. + temp_bias: float + the default value to use as temperature bias (or the starting value when + optimizing MB). Defaults to 0. + temp_bias_min: float + the minimum accepted value for the temperature bias during optimisation. + Defaults to cfg.PARAMS['temp_bias_min']. + temp_bias_max: float + the maximum accepted value for the temperature bias during optimisation. + Defaults to cfg.PARAMS['temp_bias_max']. """ # Param constraints @@ -1519,8 +1569,8 @@ def mb_calibration_from_scalar_mb(gdir, prcp_scaling_factor_min = cfg.PARAMS['prcp_scaling_factor_min'] if prcp_scaling_factor_max is None: prcp_scaling_factor_max = cfg.PARAMS['prcp_scaling_factor_max'] - if melt_f_default is None: - melt_f_default = cfg.PARAMS['melt_f_default'] + if melt_f is None: + melt_f = cfg.PARAMS['melt_f'] if temp_bias_min is None: temp_bias_min = cfg.PARAMS['temp_bias_min'] if temp_bias_max is None: @@ -1570,7 +1620,6 @@ def mb_calibration_from_scalar_mb(gdir, 'you posted!') # Ok, regardless on how we want to calibrate, we start with defaults - melt_f = melt_f_default if prcp_scaling_factor is None: if cfg.PARAMS['use_winter_prcp_factor']: # Some sanity check diff --git a/oggm/params.cfg b/oggm/params.cfg index 6cdb523b0..09d183e71 100644 --- a/oggm/params.cfg +++ b/oggm/params.cfg @@ -200,7 +200,7 @@ temp_melt = -1. # These values are given in kg m-2 K-1 day-1, for both the daily and # monthly models. We convert them back to monthly values where needed # (monthly = daily * 365 / 12) -melt_f_default = 5 +melt_f = 5 melt_f_min = 1.5 melt_f_max = 17 diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index df6bfed94..fde000a3a 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -3123,7 +3123,7 @@ def test_random(self, hef_gdir, inversion_params): dfo = hef_gdir.read_json('mb_calib') df = massbalance.mb_calibration_from_scalar_mb(hef_gdir, calibrate_param1='temp_bias', - melt_f_default=dfo['melt_f'], + melt_f=dfo['melt_f'], ref_mb=0, ref_mb_years=(1970, 2001), write_to_gdir=False) diff --git a/oggm/tests/test_prepro.py b/oggm/tests/test_prepro.py index f82c4c789..2e49f931d 100644 --- a/oggm/tests/test_prepro.py +++ b/oggm/tests/test_prepro.py @@ -1349,12 +1349,12 @@ def test_distribute_climate_historical_alps_new(self): np.testing.assert_allclose(nc_h.ref_pix_dis, nc_c.ref_pix_dis) - def test_geodetic_mb_calibration(self): + def test_mb_calibration_from_scalar_mb(self): from oggm.core.massbalance import mb_calibration_from_scalar_mb from functools import partial - mb_calibration_from_geodetic_mb = partial(mb_calibration_from_scalar_mb, - overwrite_gdir=True) + mb_calibration_from_scalar_mb = partial(mb_calibration_from_scalar_mb, + overwrite_gdir=True) hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] @@ -1372,9 +1372,9 @@ def test_geodetic_mb_calibration(self): ref_period = f'{mbdf.index[0]}-01-01_{mbdf.index[-1] + 1}-01-01' # Default is to calibrate melt_f - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period) + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period) h, w = gdir.get_inversion_flowline_hw() mb_new = massbalance.MonthlyTIModel(gdir) @@ -1388,14 +1388,14 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] == 0 - assert pdf['melt_f'] != cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] != cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] == cfg.PARAMS['prcp_scaling_factor'] # Let's calibrate on temp_bias - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='temp_bias') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='temp_bias') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['temp_mb'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1408,14 +1408,14 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] != 0 - assert pdf['melt_f'] == cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] == cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] == cfg.PARAMS['prcp_scaling_factor'] # Let's calibrate on precip - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['prcp_mb'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1428,7 +1428,7 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] == 0 - assert pdf['melt_f'] == cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] == cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] != cfg.PARAMS['prcp_scaling_factor'] # mbdf[['ref_mb', 'melt_mb', 'temp_mb', 'prcp_mb']].plot() @@ -1438,13 +1438,13 @@ def test_geodetic_mb_calibration(self): # Very positive ref_mb = 2000 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period) - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param2='temp_bias') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period) + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param2='temp_bias') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['melt_mb2'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1457,20 +1457,20 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] < 0 - assert pdf['melt_f'] != cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] != cfg.PARAMS['melt_f'] assert pdf['melt_f'] == cfg.PARAMS['melt_f_min'] assert pdf['prcp_fac'] == cfg.PARAMS['prcp_scaling_factor'] # Very negative ref_mb = -10000 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period) - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param2='temp_bias') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period) + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param2='temp_bias') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['melt_mb2'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1483,7 +1483,7 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] > 0 - assert pdf['melt_f'] != cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] != cfg.PARAMS['melt_f'] assert pdf['melt_f'] == cfg.PARAMS['melt_f_max'] assert pdf['prcp_fac'] == cfg.PARAMS['prcp_scaling_factor'] @@ -1491,15 +1491,15 @@ def test_geodetic_mb_calibration(self): # Very positive ref_mb = 3000 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac') - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['melt_mb2'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1512,21 +1512,21 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] < 0 - assert pdf['melt_f'] == cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] == cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] > cfg.PARAMS['prcp_scaling_factor'] # Very negative ref_mb = -10000 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac') - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['melt_mb2'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1539,30 +1539,30 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] > 0 - assert pdf['melt_f'] == cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] == cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] < cfg.PARAMS['prcp_scaling_factor'] # Extremely negative ref_mb = -20000 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac') with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias') - - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias', - calibrate_param3='melt_f') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias') + + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias', + calibrate_param3='melt_f') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['melt_mb3'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1575,18 +1575,18 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] == cfg.PARAMS['temp_bias_max'] - assert pdf['melt_f'] > cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] > cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] == cfg.PARAMS['prcp_scaling_factor_min'] # Unmatchable positive ref_mb = 10000 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias', - calibrate_param3='melt_f') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias', + calibrate_param3='melt_f') # Matchable positive with less range ref_mb = 1000 @@ -1595,24 +1595,24 @@ def test_geodetic_mb_calibration(self): cfg.PARAMS['prcp_scaling_factor_min'] = 2 cfg.PARAMS['prcp_scaling_factor_max'] = 3 with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac') with pytest.raises(RuntimeError): - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias') - - mb_calibration_from_geodetic_mb(gdir, - ref_mb=ref_mb, - ref_period=ref_period, - calibrate_param1='prcp_fac', - calibrate_param2='temp_bias', - calibrate_param3='melt_f') + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias') + + mb_calibration_from_scalar_mb(gdir, + ref_mb=ref_mb, + ref_period=ref_period, + calibrate_param1='prcp_fac', + calibrate_param2='temp_bias', + calibrate_param3='melt_f') mb_new = massbalance.MonthlyTIModel(gdir) mbdf['melt_mb3'] = mb_new.get_specific_mb(h, w, year=mbdf.index) @@ -1625,17 +1625,17 @@ def test_geodetic_mb_calibration(self): pdf = gdir.read_json('mb_calib') assert pdf['temp_bias'] == cfg.PARAMS['temp_bias_min'] - assert pdf['melt_f'] < cfg.PARAMS['melt_f_default'] + assert pdf['melt_f'] < cfg.PARAMS['melt_f'] assert pdf['prcp_fac'] == cfg.PARAMS['prcp_scaling_factor_max'] @pytest.mark.slow - def test_geodetic_mb_calibration_multiple_fl(self): + def test_mb_calibration_from_scalar_mb_multiple_fl(self): from oggm.core.massbalance import mb_calibration_from_scalar_mb from functools import partial - mb_calibration_from_geodetic_mb = partial(mb_calibration_from_scalar_mb, - overwrite_gdir=True) + mb_calibration_from_scalar_mb = partial(mb_calibration_from_scalar_mb, + overwrite_gdir=True) hef_file = get_demo_file('Hintereisferner_RGI5.shp') entity = gpd.read_file(hef_file).iloc[0] @@ -1654,7 +1654,7 @@ def test_geodetic_mb_calibration_multiple_fl(self): mbdf['ref_mb'] = mbdf['ANNUAL_BALANCE'] ref_mb = mbdf.ref_mb.mean() ref_period = f'{mbdf.index[0]}-01-01_{mbdf.index[-1] + 1}-01-01' - mb_calibration_from_geodetic_mb(gdir, ref_mb=ref_mb, + mb_calibration_from_scalar_mb(gdir, ref_mb=ref_mb, ref_period=ref_period) mb_new = massbalance.MonthlyTIModel(gdir) @@ -1677,7 +1677,7 @@ def test_geodetic_mb_calibration_multiple_fl(self): fls[1].surface_h -= 700 gdir.write_pickle(fls, 'inversion_flowlines') - mb_calibration_from_geodetic_mb(gdir, ref_mb=ref_mb, + mb_calibration_from_scalar_mb(gdir, ref_mb=ref_mb, ref_period=ref_period) mb_new = massbalance.MultipleFlowlineMassBalance(gdir, use_inversion_flowlines=True)