Skip to content

Commit

Permalink
Introduce scaling factor for err_dmdtda and suppress logging, both in…
Browse files Browse the repository at this point in the history 
… dynamic mu star calibration (#1488)

* added err_dmdtda_scaling_factor in dynamic mu star calibration to reduce individual uncertainty

* suppress logging of intermediate steps in dynamic mu star calibration

* switched off logging also in fallback function of dynamic mu star calibration
  • Loading branch information
@pat-schmitt
pat-schmitt committed Oct 4, 2022
1 parent c07343f commit 1a97f3b
Show file tree
Hide file tree
Showing 5 changed files with 180 additions and 36 deletions.
17 changes: 14 additions & 3 deletions oggm/cli/prepro_levels.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,8 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
add_consensus=False, start_level=None,
start_base_url=None, max_level=5, ref_tstars_base_url='',
logging_level='WORKFLOW', disable_dl_verify=False,
dynamic_spinup=False, dynamic_spinup_start_year=1979,
dynamic_spinup=False, err_dmdtda_scaling_factor=1,
dynamic_spinup_start_year=1979,
continue_on_error=True):
"""Generate the preprocessed OGGM glacier directories for this OGGM version
Expand Down Expand Up @@ -154,8 +155,11 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
a dict of parameters to override.
disable_dl_verify : bool
disable the hash verification of OGGM downloads
dynamic_spinup: str
include a dynamic spinup matching 'area' OR 'volume' at the RGI-date
dynamic_spinup : str
include a dynamic spinup matching 'area/dmdtda' OR 'volume/dmdtda' at
the RGI-date
err_dmdtda_scaling_factor : float
scaling factor to reduce individual geodetic mass balance uncertainty
dynamic_spinup_start_year : int
if dynamic_spinup is set, define the starting year for the simulation.
The default is 1979, unless the climate data starts later.
Expand Down Expand Up @@ -645,6 +649,7 @@ def _time_log():

workflow.execute_entity_task(
tasks.run_dynamic_mu_star_calibration, gdirs,
err_dmdtda_scaling_factor=err_dmdtda_scaling_factor,
ys=dynamic_spinup_start_year, ye=ye,
max_mu_star=used_max_mu_star,
kwargs_run_function={'evolution_model': evolution_model,
Expand Down Expand Up @@ -835,6 +840,11 @@ def parse_args(args):
"the RGI-date, AND mass-change from Hugonnet "
"in the period 2000-2019 (dynamic mu* "
"calibration).")
parser.add_argument('--err-dmdtda-scaling-factor', type=float, default=1,
help="scaling factor to account for correlated "
"uncertainties of geodetic mass balance "
"observations when looking at regional scale. "
"Should be smaller or equal 1.")
parser.add_argument('--dynamic-spinup-start-year', type=int, default=1979,
help="if --dynamic-spinup is set, define the starting"
"year for the simulation. The default is 1979, "
Expand Down Expand Up @@ -894,6 +904,7 @@ def parse_args(args):
ref_tstars_base_url=args.ref_tstars_base_url,
evolution_model=args.evolution_model,
dynamic_spinup=dynamic_spinup,
err_dmdtda_scaling_factor=args.err_dmdtda_scaling_factor,
dynamic_spinup_start_year=args.dynamic_spinup_start_year,
)

Expand Down
67 changes: 49 additions & 18 deletions oggm/core/dynamic_spinup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1084,17 +1084,22 @@ def dynamic_mu_star_run_with_dynamic_spinup(
define_new_mu_star_in_gdir(gdir, mu_star)

if do_inversion:
apparent_mb_from_any_mb(gdir)
# do inversion with A calibration to current volume
calibrate_inversion_from_consensus(
[gdir], apply_fs_on_mismatch=True, error_on_mismatch=False,
filter_inversion_output=True,
volume_m3_reference=local_variables['vol_m3_ref'])
with utils.DisableLogger():
apparent_mb_from_any_mb(gdir,
add_to_log_file=False, # dont write to log
)
# do inversion with A calibration to current volume
calibrate_inversion_from_consensus(
[gdir], apply_fs_on_mismatch=True, error_on_mismatch=False,
filter_inversion_output=True,
volume_m3_reference=local_variables['vol_m3_ref'],
add_to_log_file=False)

# this is used to keep the original model_flowline unchanged (-> to be able
# to conduct different dynamic calibration runs in the same gdir)
model_flowline_filesuffix = '_dyn_mu_calib'
init_present_time_glacier(gdir, filesuffix=model_flowline_filesuffix)
init_present_time_glacier(gdir, filesuffix=model_flowline_filesuffix,
add_to_log_file=False)

# Now do a dynamic spinup to match area
# do not ignore errors in dynamic spinup, so all 'bad' files are
Expand All @@ -1103,6 +1108,7 @@ def dynamic_mu_star_run_with_dynamic_spinup(
model, last_best_t_bias = run_dynamic_spinup(
gdir,
continue_on_error=False, # force to raise an error in @entity_task
add_to_log_file=False, # dont write to log file in @entity_task
init_model_fls=fls_init,
climate_input_filesuffix=climate_input_filesuffix,
evolution_model=evolution_model,
Expand Down Expand Up @@ -1136,8 +1142,9 @@ def dynamic_mu_star_run_with_dynamic_spinup(
raise RuntimeError(f'Dynamic spinup raised error! (Message: {e})')

# calculate dmdtda from previous simulation here
ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix,
path=False)
with utils.DisableLogger():
ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix,
path=False)
dmdtda_mdl = ((ds.volume.loc[yr1_ref_mb].values -
ds.volume.loc[yr0_ref_mb].values) /
gdir.rgi_area_m2 /
Expand Down Expand Up @@ -1283,11 +1290,14 @@ def dynamic_mu_star_run_with_dynamic_spinup_fallback(
if mu_star != gdir.read_json('local_mustar')['mu_star_glacierwide']:
define_new_mu_star_in_gdir(gdir, mu_star)
if do_inversion:
apparent_mb_from_any_mb(gdir)
calibrate_inversion_from_consensus(
[gdir], apply_fs_on_mismatch=True, error_on_mismatch=False,
filter_inversion_output=True,
volume_m3_reference=local_variables['vol_m3_ref'])
with utils.DisableLogger():
apparent_mb_from_any_mb(gdir,
add_to_log_file=False)
calibrate_inversion_from_consensus(
[gdir], apply_fs_on_mismatch=True, error_on_mismatch=False,
filter_inversion_output=True,
volume_m3_reference=local_variables['vol_m3_ref'],
add_to_log_file=False)
if os.path.isfile(os.path.join(gdir.dir,
'model_flowlines_dyn_mu_calib.pkl')):
os.remove(os.path.join(gdir.dir,
Expand All @@ -1309,6 +1319,7 @@ def dynamic_mu_star_run_with_dynamic_spinup_fallback(
model_end = run_dynamic_spinup(
gdir,
continue_on_error=False, # force to raise an error in @entity_task
add_to_log_file=False,
init_model_fls=fls_init,
climate_input_filesuffix=climate_input_filesuffix,
evolution_model=evolution_model,
Expand Down Expand Up @@ -1341,7 +1352,9 @@ def dynamic_mu_star_run_with_dynamic_spinup_fallback(
'try is to conduct a run until ye without a dynamic '
'spinup.')
model_end = run_from_climate_data(
gdir, min_ys=yr_clim_min, ye=ye,
gdir,
add_to_log_file=False,
min_ys=yr_clim_min, ye=ye,
output_filesuffix=output_filesuffix,
climate_input_filesuffix=climate_input_filesuffix,
store_model_geometry=store_model_geometry,
Expand Down Expand Up @@ -1438,6 +1451,7 @@ def dynamic_mu_star_run(
model = run_from_climate_data(gdir,
# force to raise an error in @entity_task
continue_on_error=False,
add_to_log_file=False,
ys=ys, ye=ye,
output_filesuffix=output_filesuffix,
init_model_fls=fls_init,
Expand All @@ -1448,8 +1462,9 @@ def dynamic_mu_star_run(
f'(Message: {e})')

# calculate dmdtda from previous simulation here
ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix,
path=False)
with utils.DisableLogger():
ds = utils.compile_run_output(gdir, input_filesuffix=output_filesuffix,
path=False)
dmdtda_mdl = ((ds.volume.loc[yr1_ref_mb].values -
ds.volume.loc[yr0_ref_mb].values) /
gdir.rgi_area_m2 /
Expand Down Expand Up @@ -1509,6 +1524,7 @@ def dynamic_mu_star_run_fallback(
model = run_from_climate_data(gdir,
# force to raise an error in @entity_task
continue_on_error=False,
add_to_log_file=False,
ys=ys, ye=ye,
output_filesuffix=output_filesuffix,
init_model_fls=fls_init,
Expand All @@ -1524,7 +1540,7 @@ def dynamic_mu_star_run_fallback(

@entity_task(log, writes=['inversion_flowlines'])
def run_dynamic_mu_star_calibration(
gdir, ref_dmdtda=None, err_ref_dmdtda=None,
gdir, ref_dmdtda=None, err_ref_dmdtda=None, err_dmdtda_scaling_factor=1,
ref_period='', ignore_hydro_months=False, min_mu_star=None,
max_mu_star=None, mu_star_max_step_length=5, maxiter=20,
ignore_errors=False, output_filesuffix='_dynamic_mu_star',
Expand Down Expand Up @@ -1563,6 +1579,20 @@ def run_dynamic_mu_star_calibration(
yr-1). Must always be a positive number. If None the data from Hugonett
2021 is used.
Default is None
err_dmdtda_scaling_factor : float
The error of the geodetic mass balance is multiplied by this factor.
When looking at more glaciers you should set this factor smaller than
1 (Default), but the smaller this factor the more glaciers will fail
during calibration. The factor is only used if ref_dmdtda = None and
err_ref_dmdtda = None.
The idea is that we reduce the uncertainty of individual observations
to count for correlated uncertainties when looking at regional or
global scales. If err_scaling_factor is 1 (Default) and you look at the
results of more than one glacier this equals that all errors are
uncorrelated. Therefore the result will be outside the uncertainty
boundaries given in Hugonett 2021 e.g. for the global estimate, because
some correlation of the individual errors is assumed during aggregation
of glaciers to regions (for more details see paper Hugonett 2021).
ref_period : str
If ref_dmdtda is None one of '2000-01-01_2010-01-01',
'2010-01-01_2020-01-01', '2000-01-01_2020-01-01'. If ref_dmdtda is
Expand Down Expand Up @@ -1697,6 +1727,7 @@ def run_dynamic_mu_star_calibration(
err_ref_dmdtda = float(df_ref_dmdtda.loc[df_ref_dmdtda['period'] ==
ref_period]['err_dmdtda'])
err_ref_dmdtda *= 1000 # kg m-2 yr-1
err_ref_dmdtda *= err_dmdtda_scaling_factor

if err_ref_dmdtda <= 0:
raise RuntimeError('The provided error for the geodetic mass-balance '
Expand Down
83 changes: 83 additions & 0 deletions oggm/tests/test_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3835,6 +3835,57 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self,
for fl_prev, fl_now in
zip(fls, gdir.read_pickle('model_flowlines'))])

# test err_dmdtda_scaling_factor (not working for volume with inversion)
if not (do_inversion and minimise_for == 'volume'):
err_dmdtda_scaling_factor = 0.2
run_dynamic_mu_star_calibration(
gdir, max_mu_star=1000.,
err_dmdtda_scaling_factor=err_dmdtda_scaling_factor,
run_function=dynamic_mu_star_run_with_dynamic_spinup,
kwargs_run_function={'minimise_for': minimise_for,
'precision_percent': precision_percent,
'precision_absolute': precision_absolute,
'do_inversion': do_inversion},
fallback_function=dynamic_mu_star_run_with_dynamic_spinup_fallback,
kwargs_fallback_function={'minimise_for': minimise_for,
'precision_percent': precision_percent,
'precision_absolute': precision_absolute,
'do_inversion': do_inversion},
output_filesuffix='_dyn_mu_calib_err_scaling',
ys=1979, ye=ye)

# check that we are matching all desired ref values
ds = utils.compile_run_output(
gdir, input_filesuffix='_dyn_mu_calib_err_scaling',
path=False)
if minimise_for == 'volume':
assert np.isclose(ds.loc[{'time': yr_rgi}].volume.values * 1e-9,
ref_value_dynamic_spinup,
rtol=precision_percent / 100,
atol=precision_absolute)
elif minimise_for == 'area':
assert np.isclose(ds.loc[{'time': yr_rgi}].area.values * 1e-6,
ref_value_dynamic_spinup,
rtol=precision_percent / 100,
atol=precision_absolute)
dmdtda_mdl_scale = ((ds.volume.loc[yr1_ref_dmdtda].values -
ds.volume.loc[yr0_ref_dmdtda].values) /
gdir.rgi_area_m2 /
(yr1_ref_dmdtda - yr0_ref_dmdtda) *
cfg.PARAMS['ice_density'])
assert np.isclose(dmdtda_mdl_scale, ref_dmdtda,
rtol=np.abs(err_ref_dmdtda *
err_dmdtda_scaling_factor / ref_dmdtda))
# check that calibration without scaling factor is outside adapted
# uncertainty (use previous result without err scaling factor), this
# tests if the scaling factor has any effect
assert not np.isclose(dmdtda_mdl, ref_dmdtda,
rtol=np.abs(err_ref_dmdtda *
err_dmdtda_scaling_factor /
ref_dmdtda))
assert gdir.get_diagnostics()['used_spinup_option'] == \
'dynamic mu_star calibration (full success)'

# test that error is raised if user provides flowlines but want to
# include inversion during dynamic mu calibration
if do_inversion:
Expand Down Expand Up @@ -4120,6 +4171,38 @@ def test_run_dynamic_mu_star_calibration_without_dynamic_spinup(self):
assert gdir.get_diagnostics()['used_spinup_option'] == \
'dynamic mu_star calibration (full success)'

# test err_dmdtda_scaling_factor
err_dmdtda_scaling_factor = 0.01
run_dynamic_mu_star_calibration(
gdir, max_mu_star=1000.,
err_dmdtda_scaling_factor=err_dmdtda_scaling_factor,
run_function=dynamic_mu_star_run,
fallback_function=dynamic_mu_star_run_fallback,
output_filesuffix='_dyn_mu_calib_err_scaling',
ys=1979, ye=ye)
# check that we are matching all desired ref values
ds = utils.compile_run_output(
gdir, input_filesuffix='_dyn_mu_calib_err_scaling',
path=False)

dmdtda_mdl_scale = ((ds.volume.loc[yr1_ref_dmdtda].values -
ds.volume.loc[yr0_ref_dmdtda].values) /
gdir.rgi_area_m2 /
(yr1_ref_dmdtda - yr0_ref_dmdtda) *
cfg.PARAMS['ice_density'])
assert np.isclose(dmdtda_mdl_scale, ref_dmdtda,
rtol=np.abs(err_ref_dmdtda *
err_dmdtda_scaling_factor / ref_dmdtda))
# check that calibration without scaling factor is outside adapted
# uncertainty (use previous result without err scaling factor), this
# tests if the scaling factor has any effect
assert not np.isclose(dmdtda_mdl, ref_dmdtda,
rtol=np.abs(err_ref_dmdtda *
err_dmdtda_scaling_factor /
ref_dmdtda))
assert gdir.get_diagnostics()['used_spinup_option'] == \
'dynamic mu_star calibration (full success)'

# tests for user provided dmdtda (always reset gdir before each test)
gdir = workflow.init_glacier_directories(
['RGI60-11.00897'], # Hintereisferner
Expand Down
3 changes: 2 additions & 1 deletion oggm/utils/_workflow.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -493,7 +493,8 @@ def _entity_task(gdir, *, reset=None, print_log=True,
if cfg.PARAMS['task_timeout'] > 0:
signal.alarm(0)
if task_name != 'gdir_to_tar':
gdir.log(task_name, task_time=ex_t)
if add_to_log_file:
gdir.log(task_name, task_time=ex_t)
except Exception as err:
# Something happened
out = None
Expand Down
Loading

0 comments on commit 1a97f3b

Please sign in to comment.