From 63918c25b08bfd1144a2622871d266d73d07eb20 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 11 Nov 2025 14:59:08 +0000 Subject: [PATCH 01/18] update bias file --- docs/_code/prepare_hef.py | 1 - oggm/cfg.py | 3 +-- oggm/cli/prepro_levels.py | 4 ---- oggm/core/massbalance.py | 10 +++------- oggm/params.cfg | 3 --- oggm/tests/funcs.py | 3 --- oggm/tests/test_benchmarks.py | 2 -- oggm/tests/test_graphics.py | 5 ----- oggm/tests/test_models.py | 3 --- oggm/tests/test_prepro.py | 5 ----- oggm/tests/test_utils.py | 6 ------ oggm/tests/test_workflow.py | 1 - oggm/utils/_downloads.py | 2 +- 13 files changed, 5 insertions(+), 43 deletions(-) diff --git a/docs/_code/prepare_hef.py b/docs/_code/prepare_hef.py index 93014ed58..0abd4e347 100644 --- a/docs/_code/prepare_hef.py +++ b/docs/_code/prepare_hef.py @@ -36,7 +36,6 @@ 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) diff --git a/oggm/cfg.py b/oggm/cfg.py index c76963a51..6eb724a58 100644 --- a/oggm/cfg.py +++ b/oggm/cfg.py @@ -532,7 +532,6 @@ def initialize_minimal(file=None, logging_level='INFO', params=None): 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)] @@ -575,7 +574,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', diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 66c543617..72ef2f938 100644 --- a/oggm/cli/prepro_levels.py +++ b/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' diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index e2f654c21..defbddbf1 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1628,11 +1628,11 @@ 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.') + raise InvalidWorkflowError('only for W5E5 data.') + bias_df = get_temp_bias_dataframe() ref_lon = climinfo['baseline_climate_ref_pix_lon'] ref_lat = climinfo['baseline_climate_ref_pix_lat'] @@ -1642,10 +1642,6 @@ def mb_calibration_from_geodetic_mb(gdir, *, temp_bias = sel_df['median_temp_bias_w_err_grouped'] 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`.') diff --git a/oggm/params.cfg b/oggm/params.cfg index eedc2e433..27377d85e 100644 --- a/oggm/params.cfg +++ b/oggm/params.cfg @@ -228,9 +228,6 @@ 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 diff --git a/oggm/tests/funcs.py b/oggm/tests/funcs.py index 5c96801fd..7639c6a56 100644 --- a/oggm/tests/funcs.py +++ b/oggm/tests/funcs.py @@ -403,7 +403,6 @@ def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None, 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 @@ -511,7 +510,6 @@ def init_columbia(reset=False): 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' @@ -552,7 +550,6 @@ def init_columbia_eb(dir_name, reset=False): 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' diff --git a/oggm/tests/test_benchmarks.py b/oggm/tests/test_benchmarks.py index 6b218f18d..ef20a1ec3 100644 --- a/oggm/tests/test_benchmarks.py +++ b/oggm/tests/test_benchmarks.py @@ -51,7 +51,6 @@ def setUp(self): 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' @@ -486,7 +485,6 @@ def setUp(self): 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' diff --git a/oggm/tests/test_graphics.py b/oggm/tests/test_graphics.py index 63818dfc2..8bbe11be1 100644 --- a/oggm/tests/test_graphics.py +++ b/oggm/tests/test_graphics.py @@ -180,7 +180,6 @@ def test_multiple_inversion(): 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 @@ -266,7 +265,6 @@ def test_multiple_models(): 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 @@ -377,7 +375,6 @@ def test_chhota_shigri(): 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') @@ -422,7 +419,6 @@ def test_ice_cap(): 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')) @@ -465,7 +461,6 @@ def test_coxe(): 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') diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index d8557a89f..da84db6a3 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -309,7 +309,6 @@ def other_glacier_cfg(): cfg.PATHS['dem_file'] = get_demo_file('srtm_oetztal.tif') cfg.PATHS['climate_file'] = get_demo_file('histalp_merged_hef.nc') 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' @@ -5341,7 +5340,6 @@ def test_hef_retreat(self, class_case_dir): cfg.PARAMS['use_multiprocessing'] = False cfg.PARAMS['min_ice_thick_for_length'] = 5 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('Hintereisferner_RGI5.shp') @@ -5438,7 +5436,6 @@ def merged_hef_cfg(class_case_dir): cfg.PARAMS['prcp_fac'] = 1.75 cfg.PARAMS['temp_melt'] = -1.75 cfg.PARAMS['use_winter_prcp_fac'] = False - cfg.PARAMS['use_temp_bias_from_file'] = False @pytest.mark.usefixtures('merged_hef_cfg') diff --git a/oggm/tests/test_prepro.py b/oggm/tests/test_prepro.py index 6668f7cf3..a58e31e0c 100644 --- a/oggm/tests/test_prepro.py +++ b/oggm/tests/test_prepro.py @@ -847,7 +847,6 @@ def setUp(self): cfg.PATHS['climate_file'] = get_demo_file('histalp_merged_hef.nc') cfg.PARAMS['baseline_climate'] = '' cfg.PARAMS['use_winter_prcp_fac'] = False - cfg.PARAMS['use_temp_bias_from_file'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -1851,7 +1850,6 @@ def setUp(self): cfg.PARAMS['baseline_climate'] = '' 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 def tearDown(self): @@ -2397,7 +2395,6 @@ def setUp(self): cfg.PATHS['working_dir'] = self.testdir cfg.PARAMS['border'] = 40 cfg.PARAMS['use_winter_prcp_fac'] = False - cfg.PARAMS['use_temp_bias_from_file'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -2516,7 +2513,6 @@ def setUp(self): cfg.PARAMS['use_multiple_flowlines'] = False cfg.PARAMS['use_tar_shapefiles'] = False cfg.PARAMS['use_winter_prcp_fac'] = False - cfg.PARAMS['use_temp_bias_from_file'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -2649,7 +2645,6 @@ def setUp(self): cfg.PATHS['climate_file'] = '' 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' diff --git a/oggm/tests/test_utils.py b/oggm/tests/test_utils.py index 5150d9cf8..dd04a668d 100644 --- a/oggm/tests/test_utils.py +++ b/oggm/tests/test_utils.py @@ -854,7 +854,6 @@ def test_start_from_level_3(self): cfg.PARAMS['prcp_fac'] = 2.5 cfg.PARAMS['use_winter_prcp_fac'] = False - cfg.PARAMS['use_temp_bias_from_file'] = False n_intersects = 0 for gdir in gdirs: @@ -1287,7 +1286,6 @@ def test_full_run_cru_centerlines(self): override_params={'geodetic_mb_period': ref_period, 'baseline_climate': 'CRU', 'use_winter_prcp_fac': False, - 'use_temp_bias_from_file': False, 'prcp_fac': 2.5, } ) @@ -1456,7 +1454,6 @@ def test_elev_bands_and_spinup_run_with_different_evolution_models(self): 'evolution_model': evolution_model, 'downstream_line_shape': downstream_line_shape, 'use_winter_prcp_fac': False, - 'use_temp_bias_from_file': False, 'prcp_fac': 2.5, }) @@ -1596,7 +1593,6 @@ def test_geodetic_per_glacier_and_massredis_run(self): params = {'geodetic_mb_period': '2000-01-01_2010-01-01', 'baseline_climate': 'CRU', 'use_winter_prcp_fac': False, - 'use_temp_bias_from_file': False, 'prcp_fac': 2.5, 'evolution_model': 'MassRedistributionCurve', 'downstream_line_shape': 'parabola', @@ -1708,7 +1704,6 @@ def test_start_from_prepro(self): override_params={'geodetic_mb_period': ref_period, 'baseline_climate': 'CRU', 'use_winter_prcp_fac': False, - 'use_temp_bias_from_file': False, 'prcp_fac': 2.5, } ) @@ -1910,7 +1905,6 @@ def test_full_benchmark_run(self): 'geodetic_mb_period': '2000-01-01_2010-01-01', 'baseline_climate': 'CRU', 'use_winter_prcp_fac': False, - 'use_temp_bias_from_file': False, 'prcp_fac': 2.5, }) diff --git a/oggm/tests/test_workflow.py b/oggm/tests/test_workflow.py index 7d0fdbadd..e41b159a6 100644 --- a/oggm/tests/test_workflow.py +++ b/oggm/tests/test_workflow.py @@ -99,7 +99,6 @@ def up_to_climate(reset=False, use_mp=None, params_file=None): cfg.PARAMS['use_kcalving_for_run'] = True cfg.PARAMS['store_model_geometry'] = True cfg.PARAMS['use_winter_prcp_fac'] = False - cfg.PARAMS['use_temp_bias_from_file'] = False cfg.PARAMS['baseline_climate'] = 'CRU' cfg.PARAMS['evolution_model'] = 'FluxBased' cfg.PARAMS['downstream_line_shape'] = 'parabola' diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 75fb00fa4..3a7c1fe95 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1340,7 +1340,7 @@ def get_temp_bias_dataframe(dataset='w5e5'): # fetch the file online base_url = ('https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' - 'w5e5_temp_bias_v2023.4.csv') + 'w5e5_rgi6_perglacier_temp_bias_v2025.6.csv') file_path = file_downloader(base_url) From f338ec56858dd7634874ce392cd989625b833de6 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 20 Nov 2025 10:59:35 +0000 Subject: [PATCH 02/18] back to old --- oggm/utils/_downloads.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 3a7c1fe95..08f7584db 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1340,7 +1340,7 @@ def get_temp_bias_dataframe(dataset='w5e5'): # fetch the file online base_url = ('https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' - 'w5e5_rgi6_perglacier_temp_bias_v2025.6.csv') + 'w5e5_rgi6_perglacier_temp_bias_v2025.6.2.csv') file_path = file_downloader(base_url) From 064ba41bbefb58f7f7965ab3d992d0663b321f77 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 2 Dec 2025 15:05:55 +0000 Subject: [PATCH 03/18] Clean up winter_precip_fac and add ERA5 --- docs/_code/prepare_hef.py | 1 - oggm/cfg.py | 7 ++--- oggm/core/massbalance.py | 51 ++++++++++++++++++----------------- oggm/params.cfg | 9 +------ oggm/tests/funcs.py | 3 --- oggm/tests/test_benchmarks.py | 2 -- oggm/tests/test_graphics.py | 5 ---- oggm/tests/test_models.py | 3 --- oggm/tests/test_prepro.py | 6 ----- oggm/tests/test_shop.py | 1 - oggm/tests/test_utils.py | 6 ----- oggm/tests/test_workflow.py | 1 - 12 files changed, 30 insertions(+), 65 deletions(-) diff --git a/docs/_code/prepare_hef.py b/docs/_code/prepare_hef.py index 0abd4e347..297b75dad 100644 --- a/docs/_code/prepare_hef.py +++ b/docs/_code/prepare_hef.py @@ -35,7 +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['prcp_fac'] = 2.5 tasks.process_histalp_data(gdir) tasks.mb_calibration_from_wgms_mb(gdir) diff --git a/oggm/cfg.py b/oggm/cfg.py index 6eb724a58..29fffc384 100644 --- a/oggm/cfg.py +++ b/oggm/cfg.py @@ -531,10 +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') - 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' @@ -585,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) diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index defbddbf1..fbe6f0217 100644 --- a/oggm/core/massbalance.py +++ b/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'], @@ -1630,9 +1642,6 @@ def mb_calibration_from_geodetic_mb(gdir, *, temp_bias = 0 if informed_threestep: climinfo = gdir.get_climate_info() - if 'w5e5' not in climinfo['baseline_climate_source'].lower(): - raise InvalidWorkflowError('only for W5E5 data.') - bias_df = get_temp_bias_dataframe() ref_lon = climinfo['baseline_climate_ref_pix_lon'] ref_lat = climinfo['baseline_climate_ref_pix_lat'] @@ -1642,9 +1651,10 @@ def mb_calibration_from_geodetic_mb(gdir, *, temp_bias = sel_df['median_temp_bias_w_err_grouped'] assert np.isfinite(temp_bias), 'Temp bias not finite?' - 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. @@ -1880,17 +1890,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 diff --git a/oggm/params.cfg b/oggm/params.cfg index 27377d85e..68ac616d9 100644 --- a/oggm/params.cfg +++ b/oggm/params.cfg @@ -230,19 +230,12 @@ temp_bias_min = -15 temp_bias_max = 15 # 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' diff --git a/oggm/tests/funcs.py b/oggm/tests/funcs.py index 7639c6a56..b20b10730 100644 --- a/oggm/tests/funcs.py +++ b/oggm/tests/funcs.py @@ -402,7 +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['evolution_model'] = 'FluxBased' cfg.PARAMS['downstream_line_shape'] = 'parabola' cfg.PARAMS['prcp_fac'] = 2.5 @@ -509,7 +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['prcp_fac'] = 2.5 cfg.PARAMS['baseline_climate'] = 'CRU' cfg.PARAMS['evolution_model'] = 'FluxBased' @@ -549,7 +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['prcp_fac'] = 2.5 cfg.PARAMS['baseline_climate'] = 'CRU' cfg.PARAMS['evolution_model'] = 'FluxBased' diff --git a/oggm/tests/test_benchmarks.py b/oggm/tests/test_benchmarks.py index ef20a1ec3..5267cc709 100644 --- a/oggm/tests/test_benchmarks.py +++ b/oggm/tests/test_benchmarks.py @@ -50,7 +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['prcp_fac'] = 2.5 cfg.PARAMS['baseline_climate'] = 'CRU' @@ -484,7 +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['prcp_fac'] = 2.5 cfg.PARAMS['baseline_climate'] = 'CRU' cfg.PARAMS['evolution_model'] = 'FluxBased' diff --git a/oggm/tests/test_graphics.py b/oggm/tests/test_graphics.py index 8bbe11be1..8c4fe95e9 100644 --- a/oggm/tests/test_graphics.py +++ b/oggm/tests/test_graphics.py @@ -179,7 +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['prcp_fac'] = 2.5 cfg.PATHS['working_dir'] = testdir @@ -264,7 +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['prcp_fac'] = 2.5 cfg.PARAMS['border'] = 40 @@ -374,7 +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['prcp_fac'] = 2.5 hef_file = get_demo_file('divides_RGI50-14.15990.shp') @@ -418,7 +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['prcp_fac'] = 2.5 df = gpd.read_file(get_demo_file('divides_RGI50-05.08389.shp')) @@ -460,7 +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['prcp_fac'] = 2.5 hef_file = get_demo_file('rgi_RGI50-01.10299.shp') diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index da84db6a3..4f33e9cd1 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -308,7 +308,6 @@ def other_glacier_cfg(): cfg.set_intersects_db(get_demo_file('rgi_intersect_oetztal.shp')) cfg.PATHS['dem_file'] = get_demo_file('srtm_oetztal.tif') cfg.PATHS['climate_file'] = get_demo_file('histalp_merged_hef.nc') - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 cfg.PARAMS['baseline_climate'] = 'CRU' @@ -5339,7 +5338,6 @@ def test_hef_retreat(self, class_case_dir): cfg.PARAMS['baseline_climate'] = '' cfg.PARAMS['use_multiprocessing'] = False cfg.PARAMS['min_ice_thick_for_length'] = 5 - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 hef_file = get_demo_file('Hintereisferner_RGI5.shp') @@ -5435,7 +5433,6 @@ def merged_hef_cfg(class_case_dir): cfg.PARAMS['border'] = 100 cfg.PARAMS['prcp_fac'] = 1.75 cfg.PARAMS['temp_melt'] = -1.75 - cfg.PARAMS['use_winter_prcp_fac'] = False @pytest.mark.usefixtures('merged_hef_cfg') diff --git a/oggm/tests/test_prepro.py b/oggm/tests/test_prepro.py index a58e31e0c..b4fb89cfe 100644 --- a/oggm/tests/test_prepro.py +++ b/oggm/tests/test_prepro.py @@ -846,7 +846,6 @@ def setUp(self): cfg.PARAMS['border'] = 10 cfg.PATHS['climate_file'] = get_demo_file('histalp_merged_hef.nc') cfg.PARAMS['baseline_climate'] = '' - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -1224,7 +1223,6 @@ def setUp(self): cfg.PARAMS['temp_bias_max'] = 10 cfg.PARAMS['prcp_fac_min'] = 0.1 cfg.PARAMS['prcp_fac_max'] = 5 - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -1849,7 +1847,6 @@ def setUp(self): cfg.PATHS['climate_file'] = get_demo_file('histalp_merged_hef.nc') cfg.PARAMS['baseline_climate'] = '' cfg.PARAMS['border'] = 10 - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -2394,7 +2391,6 @@ def setUp(self): cfg.PATHS['dem_file'] = get_demo_file('dem_RGI50-01.10299.tif') cfg.PATHS['working_dir'] = self.testdir cfg.PARAMS['border'] = 40 - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -2512,7 +2508,6 @@ def setUp(self): cfg.set_intersects_db(get_demo_file('rgi_intersect_oetztal.shp')) cfg.PARAMS['use_multiple_flowlines'] = False cfg.PARAMS['use_tar_shapefiles'] = False - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 def tearDown(self): @@ -2644,7 +2639,6 @@ def setUp(self): cfg.PATHS['dem_file'] = get_demo_file('hef_srtm.tif') cfg.PATHS['climate_file'] = '' cfg.PARAMS['border'] = 10 - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['prcp_fac'] = 2.5 cfg.PARAMS['baseline_climate'] = 'CRU' diff --git a/oggm/tests/test_shop.py b/oggm/tests/test_shop.py index d47a1ab7f..41337412b 100644 --- a/oggm/tests/test_shop.py +++ b/oggm/tests/test_shop.py @@ -445,7 +445,6 @@ def test_process_gswp3_w5e5_data(self, class_case_dir): # no gradient for GSWP3-W5E5! # test climate statistics with winter_daily_mean_prcp - # they should be computed even if cfg.PARAMS['use_winter_prcp_fac'] is False! df = utils.compile_climate_statistics([gdir], path=False, add_climate_period=[1999, 2010], add_raw_climate_statistics=True, diff --git a/oggm/tests/test_utils.py b/oggm/tests/test_utils.py index dd04a668d..31b789589 100644 --- a/oggm/tests/test_utils.py +++ b/oggm/tests/test_utils.py @@ -853,7 +853,6 @@ def test_start_from_level_3(self): prepro_border=40) cfg.PARAMS['prcp_fac'] = 2.5 - cfg.PARAMS['use_winter_prcp_fac'] = False n_intersects = 0 for gdir in gdirs: @@ -1285,7 +1284,6 @@ def test_full_run_cru_centerlines(self): centerlines=True, override_params={'geodetic_mb_period': ref_period, 'baseline_climate': 'CRU', - 'use_winter_prcp_fac': False, 'prcp_fac': 2.5, } ) @@ -1453,7 +1451,6 @@ def test_elev_bands_and_spinup_run_with_different_evolution_models(self): 'baseline_climate': 'CRU', 'evolution_model': evolution_model, 'downstream_line_shape': downstream_line_shape, - 'use_winter_prcp_fac': False, 'prcp_fac': 2.5, }) @@ -1592,7 +1589,6 @@ def test_geodetic_per_glacier_and_massredis_run(self): params = {'geodetic_mb_period': '2000-01-01_2010-01-01', 'baseline_climate': 'CRU', - 'use_winter_prcp_fac': False, 'prcp_fac': 2.5, 'evolution_model': 'MassRedistributionCurve', 'downstream_line_shape': 'parabola', @@ -1703,7 +1699,6 @@ def test_start_from_prepro(self): logging_level='INFO', max_level=5, elev_bands=True, override_params={'geodetic_mb_period': ref_period, 'baseline_climate': 'CRU', - 'use_winter_prcp_fac': False, 'prcp_fac': 2.5, } ) @@ -1904,7 +1899,6 @@ def test_full_benchmark_run(self): 'use_multiprocessing': False, 'geodetic_mb_period': '2000-01-01_2010-01-01', 'baseline_climate': 'CRU', - 'use_winter_prcp_fac': False, 'prcp_fac': 2.5, }) diff --git a/oggm/tests/test_workflow.py b/oggm/tests/test_workflow.py index e41b159a6..122d35e01 100644 --- a/oggm/tests/test_workflow.py +++ b/oggm/tests/test_workflow.py @@ -98,7 +98,6 @@ def up_to_climate(reset=False, use_mp=None, params_file=None): cfg.PARAMS['geodetic_mb_period'] = '2000-01-01_2010-01-01' cfg.PARAMS['use_kcalving_for_run'] = True cfg.PARAMS['store_model_geometry'] = True - cfg.PARAMS['use_winter_prcp_fac'] = False cfg.PARAMS['baseline_climate'] = 'CRU' cfg.PARAMS['evolution_model'] = 'FluxBased' cfg.PARAMS['downstream_line_shape'] = 'parabola' From c43e8d5817e6deb24b7db2f14aff5f82381e1d7c Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 2 Dec 2025 22:39:52 +0000 Subject: [PATCH 04/18] add ERA5 temp_bias --- oggm/cli/prepro_levels.py | 3 ++- oggm/core/massbalance.py | 10 +++++++++- oggm/shop/ecmwf.py | 1 + oggm/utils/_downloads.py | 26 ++++++++++++++++---------- 4 files changed, 28 insertions(+), 12 deletions(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 72ef2f938..0297e5439 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -629,7 +629,8 @@ def _time_log(): # Small optim to avoid concurrency utils.get_geodetic_mb_dataframe() - utils.get_temp_bias_dataframe() + utils.get_temp_bias_dataframe(dataset='w5e5') + utils.get_temp_bias_dataframe(dataset='era5') use_regional_avg = False if '_regional' in mb_calibration_strategy: diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index fbe6f0217..43e8971b7 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1642,11 +1642,19 @@ def mb_calibration_from_geodetic_mb(gdir, *, temp_bias = 0 if informed_threestep: climinfo = gdir.get_climate_info() - bias_df = get_temp_bias_dataframe() + climsource = climinfo['baseline_climate_source'] + if 'w5e5' in climsource.lower(): + bias_df = get_temp_bias_dataframe('w5e5') + elif 'era5' in climsource.lower(): + bias_df = get_temp_bias_dataframe('era5') + 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'] assert np.isfinite(temp_bias), 'Temp bias not finite?' diff --git a/oggm/shop/ecmwf.py b/oggm/shop/ecmwf.py index 0718bc4a8..5416377af 100644 --- a/oggm/shop/ecmwf.py +++ b/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: diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 08f7584db..0fef05ca6 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1315,32 +1315,38 @@ def get_geodetic_mb_dataframe(file_path=None): return df -def get_temp_bias_dataframe(dataset='w5e5'): - """Fetches the temperature bias dataframe created by the OGGM>=v16 pre-calibration - (further explained in the OGGM mass balance tutorial: - https:// tutorials.oggm.org/stable/notebooks/tutorials/massbalance_calibration.html). - The data preparation script is available at +def get_temp_bias_dataframe(dataset): + """Fetches the temperature bias dataframe. + + The dataframe was created by the OGGM>=v16 pre-calibration + (further explained in the `OGGM mass balance tutorial `_ + + The data preparation script is available at https://nbviewer.jupyter.org/urls/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/calibration/1.6.1/prepare_bias_map.ipynb The file differs between climate datasets and OGGM versions. For W5E5 and OGGM v162, it is e.g. https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/w5e5_temp_bias_v2023.4.csv - + Parameters ---------- dataset : str - climate dataset used to choose temperature bias dataframe (currently only w5e5 available) + climate dataset used to choose temperature bias dataframe + (currently only w5e5 and era5 are available) Returns ------- a DataFrame with the data. """ - if dataset != 'w5e5': + if dataset not in ['w5e5', 'era5']: raise NotImplementedError(f'No such dataset available yet: {dataset}') # fetch the file online - base_url = ('https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' - 'w5e5_rgi6_perglacier_temp_bias_v2025.6.2.csv') + base_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' + if dataset == 'w5e5': + base_url += 'w5e5_rgi6_perglacier_temp_bias_v2025.6.2.csv' + if dataset == 'era5': + base_url += 'era5_rgi6_perglacier_temp_bias_v2025.6.2.csv' file_path = file_downloader(base_url) From a4fa92199f45908b2b9e7e89d239fbf50938edc0 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 2 Dec 2025 23:05:49 +0000 Subject: [PATCH 05/18] add regional --- oggm/cli/prepro_levels.py | 1 + oggm/core/massbalance.py | 4 ++-- oggm/utils/_downloads.py | 7 ++++--- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 0297e5439..321bf6103 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -630,6 +630,7 @@ def _time_log(): # Small optim to avoid concurrency utils.get_geodetic_mb_dataframe() utils.get_temp_bias_dataframe(dataset='w5e5') + utils.get_temp_bias_dataframe(dataset='w5e5', regional=True) utils.get_temp_bias_dataframe(dataset='era5') use_regional_avg = False diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 43e8971b7..fc546860e 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1644,9 +1644,9 @@ def mb_calibration_from_geodetic_mb(gdir, *, climinfo = gdir.get_climate_info() climsource = climinfo['baseline_climate_source'] if 'w5e5' in climsource.lower(): - bias_df = get_temp_bias_dataframe('w5e5') + bias_df = get_temp_bias_dataframe('w5e5', regional=use_regional_avg) elif 'era5' in climsource.lower(): - bias_df = get_temp_bias_dataframe('era5') + bias_df = get_temp_bias_dataframe('era5', regional=use_regional_avg) else: raise InvalidWorkflowError('Dataset not suitable for ' f'informed 3-steps: {climsource}') diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 0fef05ca6..799895b7f 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1315,7 +1315,7 @@ def get_geodetic_mb_dataframe(file_path=None): return df -def get_temp_bias_dataframe(dataset): +def get_temp_bias_dataframe(dataset, regional=False): """Fetches the temperature bias dataframe. The dataframe was created by the OGGM>=v16 pre-calibration @@ -1343,10 +1343,11 @@ def get_temp_bias_dataframe(dataset): # fetch the file online base_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' + calibtype = 'regional' if regional else 'perglacier' if dataset == 'w5e5': - base_url += 'w5e5_rgi6_perglacier_temp_bias_v2025.6.2.csv' + base_url += f'w5e5_rgi6_{calibtype}_temp_bias_v2025.6.2.csv' if dataset == 'era5': - base_url += 'era5_rgi6_perglacier_temp_bias_v2025.6.2.csv' + base_url += f'era5_rgi6_{calibtype}_temp_bias_v2025.6.2.csv' file_path = file_downloader(base_url) From fdc8bfff17bfa6ab84d940eafa6cdbbabb551ae0 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Mon, 8 Dec 2025 23:22:06 +0000 Subject: [PATCH 06/18] update hugonnet --- oggm/core/massbalance.py | 9 ++++----- oggm/tests/test_utils.py | 10 ++++------ oggm/utils/_downloads.py | 16 +++++++++++----- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index fc546860e..41871a5d4 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1621,12 +1621,11 @@ 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 = get_geodetic_mb_dataframe(regional=True) 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'] + # 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] diff --git a/oggm/tests/test_utils.py b/oggm/tests/test_utils.py index 31b789589..33bae6fe9 100644 --- a/oggm/tests/test_utils.py +++ b/oggm/tests/test_utils.py @@ -1159,10 +1159,9 @@ def test_full_run_defaults(self): odf['AREA'] = df.rgi_area_km2 smb_oggm = np.average(odf['SMB'], weights=odf['AREA']) - dfh = 'table_hugonnet_regions_10yr_20yr_ar6period.csv' - dfh = pd.read_csv(utils.get_demo_file(dfh)) + dfh = utils.get_geodetic_mb_dataframe(regional=True) dfh = dfh.loc[dfh.period == '2000-01-01_2020-01-01'].set_index('reg') - smb_ref = dfh.loc[11, 'dmdtda'] + smb_ref = dfh.loc[11, 'dmdtda'] * 1000 np.testing.assert_allclose(smb_oggm, smb_ref, atol=200) # Whole Alps odf = pd.DataFrame(dfm.loc[2000:2019].mean(), columns=['SMB']) @@ -1326,10 +1325,9 @@ def test_full_run_cru_centerlines(self): odf['AREA'] = df.rgi_area_km2 smb_oggm = np.average(odf['SMB'], weights=odf['AREA']) - dfh = 'table_hugonnet_regions_10yr_20yr_ar6period.csv' - dfh = pd.read_csv(utils.get_demo_file(dfh)) + dfh = utils.get_geodetic_mb_dataframe(regional=True) dfh = dfh.loc[dfh.period == '2000-01-01_2020-01-01'].set_index('reg') - smb_ref = dfh.loc[11, 'dmdtda'] + smb_ref = dfh.loc[11, 'dmdtda'] * 1000 np.testing.assert_allclose(smb_oggm, smb_ref, atol=150) # Whole Alps odf = pd.DataFrame(dfm.loc[2000:2009].mean(), columns=['SMB']) diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 799895b7f..a6f23ef0d 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -68,7 +68,7 @@ # The given commit will be downloaded from github and used as source for # all sample data SAMPLE_DATA_GH_REPO = 'OGGM/oggm-sample-data' -SAMPLE_DATA_COMMIT = '9bfeb6dfea9513f790877819d9a6cbd2c7b61611' +SAMPLE_DATA_COMMIT = '00fca8809eeb6e087ba34ac0e3e713e7b185eca3' CHECKSUM_URL = 'https://cluster.klima.uni-bremen.de/data/downloads.sha256.hdf' CHECKSUM_VALIDATION_URL = CHECKSUM_URL + '.sha256' @@ -1270,7 +1270,7 @@ def _get_prepro_gdir_unlocked(rgi_version, rgi_id, border, prepro_level, return tar_base -def get_geodetic_mb_dataframe(file_path=None): +def get_geodetic_mb_dataframe(file_path=None, regional=False): """Fetches the reference geodetic dataframe for calibration. Currently that's the data from Hughonnet et al 2021, corrected for @@ -1282,6 +1282,8 @@ def get_geodetic_mb_dataframe(file_path=None): ---------- file_path : str in case you have your own file to parse (check the format first!) + regional : bool + to fetch the regional file instead - this is a different format! Returns ------- @@ -1291,8 +1293,12 @@ def get_geodetic_mb_dataframe(file_path=None): # fetch the file online or read custom file if file_path is None: base_url = 'https://cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb/' - file_name = 'hugonnet_2021_ds_rgi60_pergla_rates_10_20_worldwide_filled.hdf' - file_path = file_downloader(base_url + file_name) + if regional: + file_name = 'hugonnet_2021_regional_avg.csv' + file_path = file_downloader(base_url + file_name) + else: + file_name = 'hugonnet_2021_ds_rgi60_pergla_rates_10_20_worldwide_filled.hdf' + file_path = file_downloader(base_url + file_name) # Did we open it yet? if file_path in cfg.DATA: @@ -1301,7 +1307,7 @@ def get_geodetic_mb_dataframe(file_path=None): # If not let's go extension = os.path.splitext(file_path)[1] if extension == '.csv': - df = pd.read_csv(file_path, index_col=0) + df = pd.read_csv(file_path) elif extension == '.hdf': df = pd.read_hdf(file_path) From e5ba9846caff52c49c9b8750ba028219e5ff868f Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Tue, 9 Dec 2025 14:52:35 +0000 Subject: [PATCH 07/18] update bias file for regional calib --- oggm/core/massbalance.py | 7 ++++++- oggm/utils/_downloads.py | 5 +++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 41871a5d4..74f2c363d 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1655,7 +1655,12 @@ def mb_calibration_from_geodetic_mb(gdir, *, 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 cfg.PARAMS['prcp_fac'] is not None: diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index a6f23ef0d..8be7f0ff8 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1350,10 +1350,11 @@ def get_temp_bias_dataframe(dataset, regional=False): # fetch the file online base_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' calibtype = 'regional' if regional else 'perglacier' + version = '3' if regional else '2' if dataset == 'w5e5': - base_url += f'w5e5_rgi6_{calibtype}_temp_bias_v2025.6.2.csv' + base_url += f'w5e5_rgi6_{calibtype}_temp_bias_v2025.6.{version}.csv' if dataset == 'era5': - base_url += f'era5_rgi6_{calibtype}_temp_bias_v2025.6.2.csv' + base_url += f'era5_rgi6_{calibtype}_temp_bias_v2025.6.{version}.csv' file_path = file_downloader(base_url) From 3d08042c7cf0fb36bed6de9800e7c30540235a60 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Fri, 12 Dec 2025 23:26:17 +0000 Subject: [PATCH 08/18] invert params for RGI7 --- docs/api.rst | 1 + oggm/cli/prepro_levels.py | 17 ++++++++--- oggm/tests/test_prepro.py | 18 ++++++++++++ oggm/workflow.py | 60 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 92 insertions(+), 4 deletions(-) diff --git a/docs/api.rst b/docs/api.rst index 030a48191..c83dd45f5 100644 --- a/docs/api.rst +++ b/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 =============== diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 321bf6103..14229fa56 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -664,10 +664,19 @@ 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) + 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: diff --git a/oggm/tests/test_prepro.py b/oggm/tests/test_prepro.py index b4fb89cfe..44d8f3ded 100644 --- a/oggm/tests/test_prepro.py +++ b/oggm/tests/test_prepro.py @@ -2008,6 +2008,24 @@ def test_invert_hef_from_consensus(self): apply_fs_on_mismatch=True) np.testing.assert_allclose(df.vol_itmix_m3, df.vol_oggm_m3, rtol=0.01) + # Try with params + diags = gdir.get_diagnostics() + + dfo = workflow.invert_from_params(gdir, + glen_a=diags['inversion_glen_a'], + fs=diags['inversion_fs']) + np.testing.assert_allclose(df.vol_itmix_m3, dfo.vol_oggm_m3, rtol=0.01) + + dfi = pd.DataFrame() + dfi.index = [11] + dfi.loc[11, 'inversion_glen_a'] = diags['inversion_glen_a'] + dfi.loc[11, 'inversion_fs'] = diags['inversion_fs'] + + dfo = workflow.invert_from_params(gdir, params_df=dfi) + np.testing.assert_allclose(df.vol_itmix_m3, dfo.vol_oggm_m3, rtol=0.01) + + + @pytest.mark.slow def test_invert_hef_shapes(self): diff --git a/oggm/workflow.py b/oggm/workflow.py index e215a0c38..05ed0b45b 100644 --- a/oggm/workflow.py +++ b/oggm/workflow.py @@ -589,6 +589,7 @@ def calibrate_inversion_from_consensus(gdirs, ignore_missing=True, error_on_mismatch=True, filter_inversion_output=True, volume_m3_reference=None, + use_params_file=None, add_to_log_file=True): """Fit the total volume of the glaciers to the 2019 consensus estimate. @@ -721,6 +722,65 @@ def to_minimize(x): return df +@global_task(log) +def invert_from_params(gdirs, + params_df=None, + fs=None, glen_a=None, + filter_inversion_output=True, + add_to_log_file=True): + """instead of optimising the parameters, get them from a file. + + Useful e.g. for pre computed parameters for RGI7. + + Parameters + ---------- + gdirs : list of :py:class:`oggm.GlacierDirectory` objects + the glacier directories to process + params_df : str + the dataframe to use (currently regional) + glen_a : float + if params file is not provided, use this value + (defaults to cfg.params) + fs : float + if params file is not provided, use this value + filter_inversion_output : bool + whether or not to apply terminus thickness filtering on the inversion + output (needs the downstream lines to work). + + Returns + ------- + a dataframe with the individual glacier volumes + """ + + gdirs = utils.tolist(gdirs) + + df = pd.DataFrame({ + 'rgi_region': [gd.rgi_region for gd in gdirs] + }, index=[gd.rgi_id for gd in gdirs]) + df.index.name = 'rgi_id' + + if params_df is not None: + rgi_regs = set(df.rgi_region) + if len(rgi_regs) > 1: + raise InvalidParamsError('Glaciers from multiple RGI regions ' + 'are not supported.') + rgi_reg = int(rgi_regs.pop()) + glen_a = params_df.loc[rgi_reg, 'inversion_glen_a'] + fs = params_df.loc[rgi_reg, 'inversion_fs'] + + log.workflow(f'Applying A factor = {glen_a/cfg.PARAMS['glen_a']} ' + f'and fs = {fs}') + + # Compute the final volume with the correct A + inversion_tasks(gdirs, glen_a=glen_a, fs=fs, + filter_inversion_output=filter_inversion_output, + add_to_log_file=add_to_log_file) + df['vol_oggm_m3'] = execute_entity_task(tasks.get_inversion_volume, gdirs, + add_to_log_file=add_to_log_file) + return df + + + @global_task(log) def merge_glacier_tasks(gdirs, main_rgi_id=None, return_all=False, buffer=None, **kwargs): From 55104caa225c71692f5b6d595ab4838597735aaa Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Fri, 12 Dec 2025 23:29:35 +0000 Subject: [PATCH 09/18] small bug --- oggm/workflow.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/oggm/workflow.py b/oggm/workflow.py index 05ed0b45b..db879cc16 100644 --- a/oggm/workflow.py +++ b/oggm/workflow.py @@ -768,8 +768,8 @@ def invert_from_params(gdirs, glen_a = params_df.loc[rgi_reg, 'inversion_glen_a'] fs = params_df.loc[rgi_reg, 'inversion_fs'] - log.workflow(f'Applying A factor = {glen_a/cfg.PARAMS['glen_a']} ' - f'and fs = {fs}') + log.workflow(f"Applying A factor = {glen_a/cfg.PARAMS['glen_a']} " + f"and fs = {fs}") # Compute the final volume with the correct A inversion_tasks(gdirs, glen_a=glen_a, fs=fs, From 08bcb59aa25bef3b4465edd4643bb5d1ac6f82f8 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 13 Dec 2025 00:12:24 +0000 Subject: [PATCH 10/18] fix mini gdirs --- oggm/utils/_workflow.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py index da92eb9b8..6b2e5bac0 100644 --- a/oggm/utils/_workflow.py +++ b/oggm/utils/_workflow.py @@ -3970,8 +3970,8 @@ def copy_to_basedir(gdir, base_dir=None, setup='run'): New glacier directories from the copied folders """ base_dir = os.path.abspath(base_dir) - new_dir = os.path.join(base_dir, gdir.rgi_id[:8], gdir.rgi_id[:11], - gdir.rgi_id) + new_dir = os.path.join(base_dir, gdir.rgi_id[:-6], + gdir.rgi_id[:-3], gdir.rgi_id) if setup == 'run': paths = ['model_flowlines', 'inversion_params', 'outlines', 'mb_calib', 'climate_historical', 'glacier_grid', From f95c0417518846b8b95ce43b231908f5d6155885 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 13 Dec 2025 00:34:25 +0000 Subject: [PATCH 11/18] fix test --- oggm/core/massbalance.py | 7 +++++-- oggm/tests/test_benchmarks.py | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 74f2c363d..0c2fa5125 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1621,8 +1621,11 @@ def mb_calibration_from_geodetic_mb(gdir, *, # Get the reference data ref_mb_err = np.nan if use_regional_avg: - ref_mb_df = get_geodetic_mb_dataframe(regional=True) - ref_mb_df = ref_mb_df.loc[ref_mb_df.period == ref_period].set_index('reg') + 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 diff --git a/oggm/tests/test_benchmarks.py b/oggm/tests/test_benchmarks.py index 5267cc709..104452f9d 100644 --- a/oggm/tests/test_benchmarks.py +++ b/oggm/tests/test_benchmarks.py @@ -448,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) From dea80072ff4e8b7d2461b9861270bf7f56c5e0c8 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 13 Dec 2025 22:16:26 +0000 Subject: [PATCH 12/18] add RGI7 bias file --- oggm/cli/prepro_levels.py | 1 + oggm/core/massbalance.py | 8 ++++++-- oggm/utils/_downloads.py | 14 ++++++++++---- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 14229fa56..dc7b1a1a0 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -631,6 +631,7 @@ def _time_log(): utils.get_geodetic_mb_dataframe() 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='era5') use_regional_avg = False diff --git a/oggm/core/massbalance.py b/oggm/core/massbalance.py index 0c2fa5125..365705603 100644 --- a/oggm/core/massbalance.py +++ b/oggm/core/massbalance.py @@ -1646,9 +1646,13 @@ def mb_calibration_from_geodetic_mb(gdir, *, climinfo = gdir.get_climate_info() climsource = climinfo['baseline_climate_source'] if 'w5e5' in climsource.lower(): - bias_df = get_temp_bias_dataframe('w5e5', regional=use_regional_avg) + 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', regional=use_regional_avg) + 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}') diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 8be7f0ff8..f3a55f9a5 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1321,7 +1321,7 @@ def get_geodetic_mb_dataframe(file_path=None, regional=False): return df -def get_temp_bias_dataframe(dataset, regional=False): +def get_temp_bias_dataframe(dataset, regional=False, rgi_version='62'): """Fetches the temperature bias dataframe. The dataframe was created by the OGGM>=v16 pre-calibration @@ -1346,15 +1346,21 @@ def get_temp_bias_dataframe(dataset, regional=False): if dataset not in ['w5e5', 'era5']: raise NotImplementedError(f'No such dataset available yet: {dataset}') + if rgi_version not in ['62', '70G']: + raise NotImplementedError(f'RGI version not available yet: {rgi_version}') # fetch the file online base_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' calibtype = 'regional' if regional else 'perglacier' - version = '3' if regional else '2' + if rgi_version == '70G': + file_version = '1' + if rgi_version == '62': + file_version = '3' if regional else '2' + rgi_version = '6' if dataset == 'w5e5': - base_url += f'w5e5_rgi6_{calibtype}_temp_bias_v2025.6.{version}.csv' + base_url += f'w5e5_rgi{rgi_version}_{calibtype}_temp_bias_v2025.6.{file_version}.csv' if dataset == 'era5': - base_url += f'era5_rgi6_{calibtype}_temp_bias_v2025.6.{version}.csv' + base_url += f'era5_rgi{rgi_version}_{calibtype}_temp_bias_v2025.6.{file_version}.csv' file_path = file_downloader(base_url) From 2e7425dd8f4b98958c139af6d1289e932a90c6f7 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 27 Dec 2025 22:11:08 +0100 Subject: [PATCH 13/18] cap rgi date to 2019 --- oggm/utils/_workflow.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py index 6b2e5bac0..ee324cfb4 100644 --- a/oggm/utils/_workflow.py +++ b/oggm/utils/_workflow.py @@ -2925,6 +2925,10 @@ def __init__(self, rgi_entity, base_dir=None, reset=False, rgi_date = int(rgi_datestr[0:4]) if rgi_date < 0: rgi_date = RGI_DATE[self.rgi_region] + if rgi_date >= 2020: + log.warning(f'{self.rgi_id}: rgi_date {rgi_date} modified ' + 'to 2019 for workflow reasons.') + rgi_date = 2019 self.rgi_date = rgi_date # Root directory self.base_dir = os.path.normpath(base_dir) From 2186cd5118e748f2428772aabd5af2ed2604e9dc Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 27 Dec 2025 22:21:29 +0100 Subject: [PATCH 14/18] Fix test --- oggm/utils/_downloads.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index f3a55f9a5..364d31c95 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1346,6 +1346,8 @@ def get_temp_bias_dataframe(dataset, regional=False, rgi_version='62'): if dataset not in ['w5e5', 'era5']: raise NotImplementedError(f'No such dataset available yet: {dataset}') + if rgi_version == '60': + rgi_version = '62' if rgi_version not in ['62', '70G']: raise NotImplementedError(f'RGI version not available yet: {rgi_version}') From 09763695ef194ac48f4637e74bbb3f2a117cb5ba Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 1 Jan 2026 22:42:15 +0100 Subject: [PATCH 15/18] make RGI7C work --- oggm/cli/prepro_levels.py | 13 +++++ oggm/tasks.py | 1 + oggm/tests/test_prepro.py | 31 ++++++++++ oggm/workflow.py | 119 +++++++++++++++++++++++++++++++++++++- 4 files changed, 163 insertions(+), 1 deletion(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index dc7b1a1a0..2cc26fe1c 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -673,6 +673,19 @@ def _time_log(): 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'] * 1e-9 + workflow.execute_entity_task(tasks.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, diff --git a/oggm/tasks.py b/oggm/tasks.py index a261859ff..892a19fbb 100644 --- a/oggm/tasks.py +++ b/oggm/tasks.py @@ -68,3 +68,4 @@ from oggm.utils import copy_to_basedir from oggm.utils import gdir_to_tar from oggm.utils import merge_consecutive_run_outputs +from oggm.workflow import calibrate_inversion_from_volume diff --git a/oggm/tests/test_prepro.py b/oggm/tests/test_prepro.py index 44d8f3ded..5fda2ec48 100644 --- a/oggm/tests/test_prepro.py +++ b/oggm/tests/test_prepro.py @@ -2024,7 +2024,38 @@ def test_invert_hef_from_consensus(self): dfo = workflow.invert_from_params(gdir, params_df=dfi) np.testing.assert_allclose(df.vol_itmix_m3, dfo.vol_oggm_m3, rtol=0.01) + df = pd.read_hdf(utils.get_demo_file('rgi62_itmix_df.h5')) + # Works with Series + out = workflow.calibrate_inversion_from_volume(gdir, + vol_ref_m3=df['vol_itmix_m3']) + + df = df.loc[gdir.rgi_id] + np.testing.assert_allclose(df.vol_itmix_m3, + out['vol_oggm_m3'], + rtol=0.01) + assert out['fs'] == 0 + + # Works with float + out = workflow.calibrate_inversion_from_volume(gdir, + vol_ref_m3=df.vol_itmix_m3) + + np.testing.assert_allclose(df.vol_itmix_m3, + out['vol_oggm_m3'], + rtol=0.01) + assert out['fs'] == 0 + + # test user provided volume is working + delta_volume_m3 = 100000000 + user_provided_volume_m3 = df.vol_itmix_m3 - delta_volume_m3 + out = workflow.calibrate_inversion_from_volume( + gdir, + apply_fs_on_mismatch=True, + vol_ref_m3=user_provided_volume_m3) + np.testing.assert_allclose(user_provided_volume_m3, + out['vol_oggm_m3'], + rtol=0.01) + assert out['fs'] > 0 @pytest.mark.slow def test_invert_hef_shapes(self): diff --git a/oggm/workflow.py b/oggm/workflow.py index db879cc16..9a9cf4601 100644 --- a/oggm/workflow.py +++ b/oggm/workflow.py @@ -16,7 +16,7 @@ from oggm import cfg, tasks, utils from oggm.core import centerlines, flowline, climate, gis from oggm.exceptions import InvalidParamsError, InvalidWorkflowError -from oggm.utils import global_task +from oggm.utils import global_task, entity_task # MPI try: @@ -780,6 +780,123 @@ def invert_from_params(gdirs, return df +@entity_task(log, writes=['inversion_output']) +def calibrate_inversion_from_volume(gdir, vol_ref_m3=None, + fs=0, a_bounds=(0.1, 10), + apply_fs_on_mismatch=False, + error_on_mismatch=True, + filter_inversion_output=True): + """Fit the volume of a single glacier to a reference volume estimate. + + This is the entity task version of calibrate_inversion_from_consensus. + It finds the "best Glen A" to match the reference volume for a single glacier. + + Parameters + ---------- + gdir : :py:class:`oggm.GlacierDirectory` + the glacier directory to process + vol_ref_m3 : float + the reference volume in m3 to match. If float, take it, + if pd.Series, select the glacier, if None, error. + fs : float + invert with sliding (default: no) + a_bounds: tuple + factor to apply to default A + apply_fs_on_mismatch: bool + on mismatch, try to apply an arbitrary value of fs (fs = 5.7e-20 from + Oerlemans) and try to optimize A again. + error_on_mismatch: bool + sometimes the given bounds do not allow to find a zero mismatch: + this will normally raise an error, but you can switch this off, + use the closest value instead and move on. + filter_inversion_output : bool + whether or not to apply terminus thickness filtering on the inversion + output (needs the downstream lines to work). + + Returns + ------- + dict with the glacier volume and the calibrated parameters + """ + + if vol_ref_m3 is None: + raise InvalidParamsError('vol_ref_m3 must be provided (float or Series).') + + if isinstance(vol_ref_m3, pd.Series): + try: + vol_ref_m3 = vol_ref_m3.loc[gdir.rgi_id] + except KeyError: + raise InvalidParamsError(f'vol_ref_m3 series has no entry ' + f'for {gdir.rgi_id}.') + + # Optimize the diff to ref + def_a = cfg.PARAMS['inversion_glen_a'] + + if cfg.PARAMS['use_kcalving_for_inversion']: + raise NotImplementedError('Calving not implemented yet') + + def compute_vol(x): + # Run inversion tasks for this glacier + tasks.prepare_for_inversion(gdir, add_to_log_file=False) + tasks.mass_conservation_inversion(gdir, glen_a=x*def_a, fs=fs, + add_to_log_file=False) + if filter_inversion_output: + tasks.filter_inversion_output(gdir, add_to_log_file=False) + vol = tasks.get_inversion_volume(gdir, add_to_log_file=False) + return vol + + def to_minimize(x): + log.info(f'Volume calibration for {gdir.rgi_id} with ' + f'A factor: {x} and fs: {fs}') + vol = compute_vol(x) + return vol_ref_m3 - vol + + try: + out_fac, r = optimization.brentq(to_minimize, *a_bounds, + rtol=1e-2, + full_output=True) + if r.converged: + log.info(f'calibrate_inversion_from_volume for {gdir.rgi_id} ' + f'converged after {r.iterations} iterations and fs={fs}. The ' + f'resulting Glen A factor is {out_fac}.') + else: + raise ValueError('Unexpected error in optimization.brentq') + except ValueError: + # Ok can't find an A. Log for debug: + vol1 = compute_vol(a_bounds[0]) + vol2 = compute_vol(a_bounds[1]) + msg = (f'calibration from volume CAN\'T converge for {gdir.rgi_id} with fs={fs}.\n' + f'Bound values (m3):\nRef={vol_ref_m3:.0f} OGGM={vol1:.0f} for A factor {a_bounds[0]}\n' + f'Ref={vol_ref_m3:.0f} OGGM={vol2:.0f} for A factor {a_bounds[1]}') + if apply_fs_on_mismatch and fs == 0 and vol2 > vol_ref_m3: + return calibrate_inversion_from_volume(gdir, + vol_ref_m3=vol_ref_m3, + fs=5.7e-20, a_bounds=a_bounds, + apply_fs_on_mismatch=False, + error_on_mismatch=error_on_mismatch, + filter_inversion_output=filter_inversion_output) + if error_on_mismatch: + raise ValueError(msg) + + out_fac = a_bounds[int(abs(vol_ref_m3 - vol1) > + abs(vol_ref_m3 - vol2))] + log.info(msg) + log.info(f'We use A factor = {out_fac} and fs = {fs} and move on.') + + # Compute the final volume with the correct A + tasks.prepare_for_inversion(gdir) + tasks.mass_conservation_inversion(gdir, glen_a=out_fac*def_a, fs=fs) + if filter_inversion_output: + tasks.filter_inversion_output(gdir) + + final_vol = tasks.get_inversion_volume(gdir) + + return { + 'vol_oggm_m3': final_vol, + 'glen_a': out_fac * def_a, + 'fs': fs, + 'a_factor': out_fac + } + @global_task(log) def merge_glacier_tasks(gdirs, main_rgi_id=None, return_all=False, buffer=None, From 1c65cff9bb87b60313b1cc637a7303d99b1e6499 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Thu, 1 Jan 2026 22:47:50 +0100 Subject: [PATCH 16/18] circular import --- oggm/cli/prepro_levels.py | 2 +- oggm/tasks.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index 2cc26fe1c..dcc5799c3 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -680,7 +680,7 @@ def _time_log(): # 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'] * 1e-9 - workflow.execute_entity_task(tasks.calibrate_inversion_from_volume, + workflow.execute_entity_task(workflow.calibrate_inversion_from_volume, gdirs, vol_ref_m3=ref_vol_df, apply_fs_on_mismatch=True, diff --git a/oggm/tasks.py b/oggm/tasks.py index 892a19fbb..a261859ff 100644 --- a/oggm/tasks.py +++ b/oggm/tasks.py @@ -68,4 +68,3 @@ from oggm.utils import copy_to_basedir from oggm.utils import gdir_to_tar from oggm.utils import merge_consecutive_run_outputs -from oggm.workflow import calibrate_inversion_from_volume From e283453127dda88196b19ad18156bcffdc13c7fb Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 3 Jan 2026 23:13:59 +0100 Subject: [PATCH 17/18] add RGI7C year --- oggm/cli/prepro_levels.py | 2 ++ oggm/utils/_downloads.py | 4 ++-- oggm/utils/_workflow.py | 19 ++++++++++++++++++- 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index dcc5799c3..e47daccf8 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -629,9 +629,11 @@ def _time_log(): # Small optim to avoid concurrency utils.get_geodetic_mb_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 diff --git a/oggm/utils/_downloads.py b/oggm/utils/_downloads.py index 364d31c95..d07358a0e 100644 --- a/oggm/utils/_downloads.py +++ b/oggm/utils/_downloads.py @@ -1348,13 +1348,13 @@ def get_temp_bias_dataframe(dataset, regional=False, rgi_version='62'): raise NotImplementedError(f'No such dataset available yet: {dataset}') if rgi_version == '60': rgi_version = '62' - if rgi_version not in ['62', '70G']: + if rgi_version not in ['62', '70G', '70C']: raise NotImplementedError(f'RGI version not available yet: {rgi_version}') # fetch the file online base_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.6/' calibtype = 'regional' if regional else 'perglacier' - if rgi_version == '70G': + if rgi_version in ['70G', '70C']: file_version = '1' if rgi_version == '62': file_version = '3' if regional else '2' diff --git a/oggm/utils/_workflow.py b/oggm/utils/_workflow.py index ee324cfb4..0efbec3df 100644 --- a/oggm/utils/_workflow.py +++ b/oggm/utils/_workflow.py @@ -630,6 +630,22 @@ def get_ref_mb_glaciers(gdirs, y0=None, y1=None): return ref_gdirs +def get_rgi70C_year(rgi_id): + """Temporary function to fetch the rgi outline year for RGI70C ids. + """ + + key = 'RGI70C_rgi_year' + if key not in cfg.DATA: + from oggm.utils._downloads import file_downloader + fp = file_downloader('https://cluster.klima.uni-bremen.de/~oggm/' + 'ref_mb_params/oggm_v1.6/inv_rgi7/' + 'rgi7c_rgi_year_2025.1.csv') + + cfg.DATA[key] = pd.read_csv(fp, index_col=0)['rgi_year'] + + return int(cfg.DATA[key].loc[rgi_id]) + + def _chaikins_corner_cutting(line, refinements=5): """Some magic here. @@ -2788,7 +2804,8 @@ def __init__(self, rgi_entity, base_dir=None, reset=False, if is_glacier_complex: rgi_entity['glac_name'] = '' - rgi_entity['src_date'] = '2000-01-01 00:00:00' + rgi_year = get_rgi70C_year(self.rgi_id) + rgi_entity['src_date'] = f'{rgi_year}-01-01 00:00:00' if 'dem_source' not in rgi_entity: rgi_entity['dem_source'] = None rgi_entity['term_type'] = 9 From 9a5b7cecf9e24aa097007cc5c00a6a83b6f6b238 Mon Sep 17 00:00:00 2001 From: Fabien Maussion Date: Sat, 3 Jan 2026 23:29:29 +0100 Subject: [PATCH 18/18] order of magnitude --- oggm/cli/prepro_levels.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/oggm/cli/prepro_levels.py b/oggm/cli/prepro_levels.py index e47daccf8..0e390cadf 100644 --- a/oggm/cli/prepro_levels.py +++ b/oggm/cli/prepro_levels.py @@ -681,7 +681,7 @@ def _time_log(): 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'] * 1e-9 + 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,