update bias file

docs/_code/prepare_hef.py

@@ -35,8 +35,6 @@
 tasks.compute_downstream_line(gdir)
 tasks.compute_downstream_bedshape(gdir)
 cfg.PARAMS['baseline_climate'] = 'HISTALP'
-cfg.PARAMS['use_winter_prcp_fac'] = False
-cfg.PARAMS['use_temp_bias_from_file'] = False
 cfg.PARAMS['prcp_fac'] = 2.5
 tasks.process_histalp_data(gdir)
 tasks.mb_calibration_from_wgms_mb(gdir)

docs/api.rst

@@ -32,6 +32,7 @@ Tools to set-up and run OGGM.
     workflow.inversion_tasks
     workflow.merge_glacier_tasks
     workflow.calibrate_inversion_from_consensus
+    workflow.invert_from_params
 
 Troubleshooting
 ===============

oggm/cfg.py

@@ -531,11 +531,7 @@ def initialize_minimal(file=None, logging_level='INFO', params=None):
     PARAMS['hydro_month_nh'] = cp.as_int('hydro_month_nh')
     PARAMS['hydro_month_sh'] = cp.as_int('hydro_month_sh')
     PARAMS['geodetic_mb_period'] = cp['geodetic_mb_period']
-    PARAMS['use_winter_prcp_fac'] = cp.as_bool('use_winter_prcp_fac')
-    PARAMS['use_temp_bias_from_file'] = cp.as_bool('use_temp_bias_from_file')
 
-    k = 'winter_prcp_fac_ab'
-    PARAMS[k] = [float(vk) for vk in cp.as_list(k)]
     k = 'ref_mb_valid_window'
     PARAMS[k] = [int(vk) for vk in cp.as_list(k)]
     k = 'free_board_marine_terminating'
@@ -575,7 +571,7 @@ def initialize_minimal(file=None, logging_level='INFO', params=None):
            'grid_dx_method', 'compress_climate_netcdf', 'by_bin_dx',
            'mp_processes', 'use_multiprocessing', 'clip_dem_to_zero',
            'topo_interp', 'use_compression', 'bed_shape', 'continue_on_error',
-           'use_multiple_flowlines', 'border', 'use_temp_bias_from_file',
+           'use_multiple_flowlines', 'border',
            'mpi_recv_buf_size', 'map_proj', 'evolution_model',
            'hydro_month_sh', 'hydro_month_nh', 'by_bin_bins',
            'use_intersects', 'filter_min_slope', 'clip_tidewater_border',
@@ -586,9 +582,9 @@ def initialize_minimal(file=None, logging_level='INFO', params=None):
            'free_board_marine_terminating', 'use_kcalving_for_inversion',
            'error_when_glacier_reaches_boundaries', 'glacier_length_method',
            'use_inversion_params_for_run',
-           'tidewater_type', 'store_model_geometry', 'use_winter_prcp_fac',
+           'tidewater_type', 'store_model_geometry',
            'store_diagnostic_variables', 'store_fl_diagnostic_variables',
-           'geodetic_mb_period', 'store_fl_diagnostics', 'winter_prcp_fac_ab',
+           'geodetic_mb_period', 'store_fl_diagnostics',
            'prcp_fac', 'downstream_line_shape', 'keep_multipolygon_outlines']
     for k in ltr:
         cp.pop(k, None)

oggm/cli/prepro_levels.py

@@ -245,10 +245,6 @@ def _time_log():
     # Set to True for operational runs
     override_params['continue_on_error'] = continue_on_error
 
-    # Do not use bias file if user wants melt_temp only
-    if mb_calibration_strategy in ['melt_temp', 'temp_melt']:
-        override_params['use_temp_bias_from_file'] = False
-
     # For centerlines we have to change the default evolution model and bed
     if centerlines:
         override_params['downstream_line_shape'] = 'parabola'
@@ -633,7 +629,12 @@ def _time_log():
 
         # Small optim to avoid concurrency
         utils.get_geodetic_mb_dataframe()
-        utils.get_temp_bias_dataframe()
+        utils.get_rgi70C_year('RGI2000-v7.0-C-01-00001')
+        utils.get_temp_bias_dataframe(dataset='w5e5')
+        utils.get_temp_bias_dataframe(dataset='w5e5', regional=True)
+        utils.get_temp_bias_dataframe(dataset='w5e5', rgi_version='70G', regional=True)
+        utils.get_temp_bias_dataframe(dataset='w5e5', rgi_version='70C', regional=True)
+        utils.get_temp_bias_dataframe(dataset='era5')
 
         use_regional_avg = False
         if '_regional' in mb_calibration_strategy:
@@ -666,10 +667,32 @@ def _time_log():
 
             # Inversion: we match the consensus
             filter = border >= 20
-            workflow.calibrate_inversion_from_consensus(gdirs,
-                                                        apply_fs_on_mismatch=True,
-                                                        error_on_mismatch=False,
-                                                        filter_inversion_output=filter)
+
+            if rgi_version == '70G':
+                purl = ('https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/'
+                        'oggm_v1.6/inv_rgi7/rgi6_regional_inv_params_2025.1.csv')
+                params_df = pd.read_csv(utils.file_downloader(purl), index_col=0)
+                workflow.invert_from_params(gdirs,
+                                            params_df=params_df,
+                                            filter_inversion_output=filter)
+            elif rgi_version == '70C':
+                purl = ('https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/'
+                        'oggm_v1.6/inv_rgi7/rgi7c_glacier_inv_ref_2025.1.csv')
+                ref_vol_df = pd.read_csv(utils.file_downloader(purl), index_col=0)
+                # Small optim
+                ref_vol_df = ref_vol_df.loc[ref_vol_df.rgi_region == int(rgi_reg)]
+                ref_vol_df = ref_vol_df['inv_volume_km3'] * 1e9
+                workflow.execute_entity_task(workflow.calibrate_inversion_from_volume,
+                                             gdirs,
+                                             vol_ref_m3=ref_vol_df,
+                                             apply_fs_on_mismatch=True,
+                                             error_on_mismatch=False,
+                                             filter_inversion_output=filter)
+            else:
+                workflow.calibrate_inversion_from_consensus(gdirs,
+                                                            apply_fs_on_mismatch=True,
+                                                            error_on_mismatch=False,
+                                                            filter_inversion_output=filter)
 
             # We get ready for modelling
             if border >= 20:

oggm/core/massbalance.py

@@ -1450,13 +1450,11 @@ def calving_mb(gdir):
 
 
 def decide_winter_precip_factor(gdir):
-    """Utility function to decide on a precip factor based on winter precip."""
+    """Utility function to decide on a precip factor based on winter precip.
 
-    # We have to decide on a precip factor
-    if 'W5E5' not in cfg.PARAMS['baseline_climate']:
-        raise InvalidWorkflowError('prcp_fac from_winter_prcp is only '
-                                   'compatible with the W5E5 climate '
-                                   'dataset!')
+    The values here are hardcoded as OGGM evolves - there should be an
+    easy way to change it if people need more flexibility one day.
+    """
 
     # get non-corrected winter daily mean prcp (kg m-2 day-1)
     # it is easier to get this directly from the raw climate files
@@ -1480,10 +1478,24 @@ def decide_winter_precip_factor(gdir):
         assert len(ds_pr_winter.time) == 41 * 7, text
         w_prcp = float((ds_pr_winter / ds_pr_winter.time.dt.daysinmonth).mean())
 
-    # from MB sandbox calibration to winter MB
-    # using t_melt=-1, cte lapse rate, monthly resolution
-    a, b = cfg.PARAMS['winter_prcp_fac_ab']
-    prcp_fac = a * np.log(w_prcp) + b
+    climsource = gdir.get_climate_info()['baseline_climate_source']
+    if 'w5e5' in climsource.lower():
+        # from OGGM calibration to winter MB, LOG
+        # repeated by Lily in November 2025 with newest gdirs
+        # using t_melt=-1, cte lapse rate, monthly resolution
+        a, b = -1.0614, 3.9200
+        prcp_fac = a * np.log(w_prcp) + b
+    elif 'era5' in climsource.lower():
+        # from OGGM calibration to winter MB, LINEAR
+        # repeated by Lily in November 2025 with newest gdirs
+        # using t_melt=-1, cte lapse rate, monthly resolution
+        a, b = -0.09078476, 2.43505368
+        prcp_fac = a * w_prcp + b
+    else:
+        msg = (f'Baseline climate {climsource} not suitable for'
+               'decide_winter_precip_factor(). Set prcp_fac.')
+        raise InvalidWorkflowError(msg)
+
     # don't allow extremely low/high prcp. factors!!!
     return clip_scalar(prcp_fac,
                        cfg.PARAMS['prcp_fac_min'],
@@ -1609,12 +1621,14 @@ def mb_calibration_from_geodetic_mb(gdir, *,
     # Get the reference data
     ref_mb_err = np.nan
     if use_regional_avg:
-        ref_mb_df = 'table_hugonnet_regions_10yr_20yr_ar6period.csv'
-        ref_mb_df = pd.read_csv(get_demo_file(ref_mb_df))
-        ref_mb_df = ref_mb_df.loc[ref_mb_df.period == ref_period].set_index('reg')
-        # dmdtda already in kg m-2 yr-1
-        ref_mb = ref_mb_df.loc[int(gdir.rgi_region), 'dmdtda']
-        ref_mb_err = ref_mb_df.loc[int(gdir.rgi_region), 'err_dmdtda']
+        ref_mb_df_o = get_geodetic_mb_dataframe(regional=True)
+        ref_mb_df = ref_mb_df_o.loc[ref_mb_df_o.period == ref_period].set_index('reg')
+        if len(ref_mb_df) == 0:
+            raise InvalidParamsError(f'Ref period {ref_period} not found in file: '
+                                     f'{ref_mb_df_o.period.unique()}')
+        # dmdtda: in meters water-equivalent per year -> we convert to kg m-2 yr-1
+        ref_mb = ref_mb_df.loc[int(gdir.rgi_region), 'dmdtda'] * 1000
+        ref_mb_err = ref_mb_df.loc[int(gdir.rgi_region), 'err_dmdtda'] * 1000
     else:
         try:
             ref_mb_df = get_geodetic_mb_dataframe().loc[gdir.rgi_id]
@@ -1628,27 +1642,38 @@ def mb_calibration_from_geodetic_mb(gdir, *,
             ref_mb = override_missing
 
     temp_bias = 0
-    if cfg.PARAMS['use_temp_bias_from_file']:
+    if informed_threestep:
         climinfo = gdir.get_climate_info()
-        if 'w5e5' not in climinfo['baseline_climate_source'].lower():
-            raise InvalidWorkflowError('use_temp_bias_from_file currently '
-                                       'only available for W5E5 data.')
-        bias_df = get_temp_bias_dataframe()
+        climsource = climinfo['baseline_climate_source']
+        if 'w5e5' in climsource.lower():
+            bias_df = get_temp_bias_dataframe('w5e5',
+                                              rgi_version=gdir.rgi_version,
+                                              regional=use_regional_avg)
+        elif 'era5' in climsource.lower():
+            bias_df = get_temp_bias_dataframe('era5',
+                                              rgi_version=gdir.rgi_version,
+                                              regional=use_regional_avg)
+        else:
+            raise InvalidWorkflowError('Dataset not suitable for '
+                                       f'informed 3-steps: {climsource}')
         ref_lon = climinfo['baseline_climate_ref_pix_lon']
         ref_lat = climinfo['baseline_climate_ref_pix_lat']
         # Take nearest
         dis = ((bias_df.lon_val - ref_lon)**2 + (bias_df.lat_val - ref_lat)**2)**0.5
+        assert dis.min() < 1, 'Somethings wrong with lons'
         sel_df = bias_df.iloc[np.argmin(dis)]
-        temp_bias = sel_df['median_temp_bias_w_err_grouped']
+        # Which bias central value to use?
+        if use_regional_avg:
+            centralval = 'median_temp_bias_w_area_grouped'
+        else:
+            centralval = 'median_temp_bias_w_err_grouped'
+        temp_bias = sel_df[centralval]
         assert np.isfinite(temp_bias), 'Temp bias not finite?'
 
-    if informed_threestep:
-        if not cfg.PARAMS['use_temp_bias_from_file']:
-            raise InvalidParamsError('With `informed_threestep` you need to '
-                                     'set `use_temp_bias_from_file`.')
-        if not cfg.PARAMS['use_winter_prcp_fac']:
-            raise InvalidParamsError('With `informed_threestep` you need to '
-                                     'set `use_winter_prcp_fac`.')
+        if cfg.PARAMS['prcp_fac'] is not None:
+            raise InvalidParamsError('With `informed_threestep` you cannot use '
+                                     'a preset prcp_fac - we need to rely on '
+                                     'decide_winter_precip_factor().')
 
         # Some magic heuristics - we just decide to calibrate
         # precip -> melt_f -> temp but informed by previous data.
@@ -1884,17 +1909,10 @@ def mb_calibration_from_scalar_mb(gdir, *,
         melt_f = cfg.PARAMS['melt_f']
 
     if prcp_fac is None:
-        if cfg.PARAMS['use_winter_prcp_fac']:
-            # Some sanity check
-            if cfg.PARAMS['prcp_fac'] is not None:
-                raise InvalidWorkflowError("Set PARAMS['prcp_fac'] to None "
-                                           "if using PARAMS['winter_prcp_factor'].")
+        if cfg.PARAMS['prcp_fac'] is None:
             prcp_fac = decide_winter_precip_factor(gdir)
         else:
             prcp_fac = cfg.PARAMS['prcp_fac']
-            if prcp_fac is None:
-                raise InvalidWorkflowError("Set either PARAMS['use_winter_prcp_fac'] "
-                                           "or PARAMS['winter_prcp_factor'].")
 
     if temp_bias is None:
         temp_bias = 0

oggm/params.cfg

@@ -228,24 +228,14 @@ melt_f_max = 17
 # The default range is very large and should be used with caution
 temp_bias_min = -15
 temp_bias_max = 15
-# For W5E5 and calibration on geodetic data from Hugonnet, use the precomputed
-# bias values based on grid point analysis (see documentation)
-use_temp_bias_from_file = True
 
 # precipitation correction: set to a float for a constant scaling factor
-# Needs to be set empty if use_winter_prcp_fac is True
+# Set to None to rely on decide_winter_precip_factor() (see tutorials)
 prcp_fac =
 # For calibration - set a range
 prcp_fac_min = 0.1
 prcp_fac_max = 10
 
-# Use a precipitation dependent factor (unique per glacier)
-# The values below have been calibrated on W5E5 data and currently
-# we enforce using it with these parameters. Note that the
-# prcp range above is used to constrain this
-use_winter_prcp_fac = True
-winter_prcp_fac_ab = -1.0614, 3.9200
-
 # When matching geodetic MB on a glacier per glacier basis, which period
 # to use. Available in the current data are:
 # '2000-01-01_2010-01-01', '2010-01-01_2020-01-01', '2000-01-01_2020-01-01'

oggm/shop/ecmwf.py

@@ -158,6 +158,7 @@ def process_ecmwf_data(gdir, dataset=None, ensemble_member=0,
         temp = ds['t2m'].data - 273.15
         time = ds.time.data
         ref_lon = float(ds['longitude'])
+        # That's actually not how it should be - but it's too late to change
         ref_lon = ref_lon - 360 if ref_lon > 180 else ref_lon
         ref_lat = float(ds['latitude'])
     with xr.open_dataset(get_ecmwf_file(dataset, 'pre')) as ds:

oggm/tests/funcs.py

@@ -402,8 +402,6 @@ def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None,
     cfg.PATHS['working_dir'] = testdir
     cfg.PARAMS['trapezoid_lambdas'] = 1
     cfg.PARAMS['border'] = border
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['evolution_model'] = 'FluxBased'
     cfg.PARAMS['downstream_line_shape'] = 'parabola'
     cfg.PARAMS['prcp_fac'] = 2.5
@@ -510,8 +508,6 @@ def init_columbia(reset=False):
     cfg.PARAMS['border'] = 10
     cfg.PARAMS['use_kcalving_for_inversion'] = True
     cfg.PARAMS['use_kcalving_for_run'] = True
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['baseline_climate'] = 'CRU'
     cfg.PARAMS['evolution_model'] = 'FluxBased'
@@ -551,8 +547,6 @@ def init_columbia_eb(dir_name, reset=False):
     cfg.PARAMS['border'] = 10
     cfg.PARAMS['use_kcalving_for_inversion'] = True
     cfg.PARAMS['use_kcalving_for_run'] = True
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['baseline_climate'] = 'CRU'
     cfg.PARAMS['evolution_model'] = 'FluxBased'

oggm/tests/test_benchmarks.py

@@ -50,8 +50,6 @@ def setUp(self):
         cfg.PATHS['working_dir'] = self.testdir
         cfg.PATHS['dem_file'] = get_demo_file('dem_SouthGlacier.tif')
         cfg.PARAMS['border'] = 10
-        cfg.PARAMS['use_winter_prcp_fac'] = False
-        cfg.PARAMS['use_temp_bias_from_file'] = False
         cfg.PARAMS['prcp_fac'] = 2.5
         cfg.PARAMS['baseline_climate'] = 'CRU'
 
@@ -450,7 +448,7 @@ def test_workflow(self):
         execute_entity_task(tasks.mb_calibration_from_geodetic_mb, gdirs,
                             use_regional_avg=True,
                             overwrite_gdir=True,
-                            ref_period='2000-01-01_2010-01-01')
+                            ref_period='2000-01-01_2005-01-01')
         df['region'] = utils.compile_fixed_geometry_mass_balance(gdirs)['RGI60-01.16195']
         assert 0.99 < df.corr().iloc[0, 1] < 1
         assert np.all(df.std() > 450)
@@ -485,8 +483,6 @@ def setUp(self):
         cfg.PATHS['working_dir'] = self.testdir
         cfg.PARAMS['use_kcalving_for_inversion'] = True
         cfg.PARAMS['use_kcalving_for_run'] = True
-        cfg.PARAMS['use_winter_prcp_fac'] = False
-        cfg.PARAMS['use_temp_bias_from_file'] = False
         cfg.PARAMS['prcp_fac'] = 2.5
         cfg.PARAMS['baseline_climate'] = 'CRU'
         cfg.PARAMS['evolution_model'] = 'FluxBased'

oggm/tests/test_graphics.py

@@ -179,8 +179,6 @@ def test_multiple_inversion():
     cfg.PARAMS['border'] = 40
     cfg.PARAMS['baseline_climate'] = 'CUSTOM'
     cfg.PARAMS['trapezoid_lambdas'] = 1
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PATHS['working_dir'] = testdir
 
@@ -265,8 +263,6 @@ def test_multiple_models():
     cfg.PATHS['working_dir'] = testdir
     cfg.PARAMS['baseline_climate'] = 'CUSTOM'
     cfg.PARAMS['trapezoid_lambdas'] = 1
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
     cfg.PARAMS['border'] = 40
 
@@ -376,8 +372,6 @@ def test_chhota_shigri():
     cfg.PARAMS['use_intersects'] = False
     cfg.PATHS['working_dir'] = testdir
     cfg.PARAMS['trapezoid_lambdas'] = 1
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
 
     hef_file = get_demo_file('divides_RGI50-14.15990.shp')
@@ -421,8 +415,6 @@ def test_ice_cap():
     cfg.PARAMS['border'] = 60
     cfg.PATHS['working_dir'] = testdir
     cfg.PARAMS['trapezoid_lambdas'] = 1
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
 
     df = gpd.read_file(get_demo_file('divides_RGI50-05.08389.shp'))
@@ -464,8 +456,6 @@ def test_coxe():
     cfg.PARAMS['use_kcalving_for_inversion'] = True
     cfg.PARAMS['use_kcalving_for_run'] = True
     cfg.PARAMS['trapezoid_lambdas'] = 1
-    cfg.PARAMS['use_winter_prcp_fac'] = False
-    cfg.PARAMS['use_temp_bias_from_file'] = False
     cfg.PARAMS['prcp_fac'] = 2.5
 
     hef_file = get_demo_file('rgi_RGI50-01.10299.shp')