Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Some bug fixes for dynamic mu star calibration #1451

Merged
merged 6 commits into from
Aug 24, 2022
158 changes: 107 additions & 51 deletions oggm/core/dynamic_spinup.py

Large diffs are not rendered by default.

14 changes: 14 additions & 0 deletions oggm/core/flowline.py
Original file line number Diff line number Diff line change
Expand Up @@ -933,6 +933,7 @@ def run_until_and_store(self, y1,
stop_criterion=None,
fixed_geometry_spinup_yr=None,
dynamic_spinup_min_ice_thick=None,
make_compatible=False
):
"""Runs the model and returns intermediate steps in xarray datasets.

Expand Down Expand Up @@ -991,6 +992,13 @@ def run_until_and_store(self, y1,
area or the total volume. This is useful to smooth out yearly
fluctuations when matching to observations. The names of this new
variables include the suffix _min_h (e.g. 'area_m2_min_h')
make_compatible : bool
if set to true this will add all variables to the resulting dataset
so it could be combined with any other one. This is necessary if
different spinup methods are used. For example if using the dynamic
spinup and setting fixed geoemtry spinup as fallback, the variable
'is_fixed_geometry_spinup' must be added to the dynamic spinup so
it is possible to compile both glaciers together.
Returns
-------
geom_ds : xarray.Dataset or None
Expand Down Expand Up @@ -1161,6 +1169,12 @@ def run_until_and_store(self, y1,
desc = 'Part of the series which are spinup'
diag_ds['is_fixed_geometry_spinup'].attrs['description'] = desc
diag_ds['is_fixed_geometry_spinup'].attrs['unit'] = '-'
elif make_compatible:
is_spinup_time = np.full(len(monthly_time), False, dtype=np.bool)
diag_ds['is_fixed_geometry_spinup'] = ('time', is_spinup_time)
desc = 'Part of the series which are spinup'
diag_ds['is_fixed_geometry_spinup'].attrs['description'] = desc
diag_ds['is_fixed_geometry_spinup'].attrs['unit'] = '-'

fl_diag_dss = None
if do_fl_diag:
Expand Down
56 changes: 28 additions & 28 deletions oggm/tests/test_models.py
Original file line number Diff line number Diff line change
Expand Up @@ -3774,13 +3774,13 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self,
gdir, max_mu_star=1000.,
run_function=dynamic_mu_star_run_with_dynamic_spinup,
kwargs_run_function={'minimise_for': minimise_for,
'precision_percent_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'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_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'precision_percent': precision_percent,
'precision_absolute': precision_absolute,
'do_inversion': do_inversion},
output_filesuffix='_dyn_mu_calib_spinup_inversion',
ys=1979, ye=ye)
Expand Down Expand Up @@ -3834,13 +3834,13 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self,
gdir, max_mu_star=1000.,
run_function=dynamic_mu_star_run_with_dynamic_spinup,
kwargs_run_function={'minimise_for': minimise_for,
'precision_percent_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'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_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'precision_percent': precision_percent,
'precision_absolute': precision_absolute,
'do_inversion': do_inversion},
output_filesuffix='_dyn_mu_calib_spinup_inversion',
ys=1979, ye=ye, init_model_fls=fls)
Expand All @@ -3860,13 +3860,13 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self,
err_ref_dmdtda=err_ref_dmdtda + delta_err_ref_dmdtda,
run_function=dynamic_mu_star_run_with_dynamic_spinup,
kwargs_run_function={'minimise_for': minimise_for,
'precision_percent_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'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_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'precision_percent': precision_percent,
'precision_absolute': precision_absolute,
'do_inversion': do_inversion},
output_filesuffix='_dyn_mu_calib_spinup_inversion_user_dmdtda',
ys=1979, ye=ye)
Expand Down Expand Up @@ -3902,25 +3902,25 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self,
output_filesuffix='_dyn_mu_calib_spinup_inversion_error',
ignore_errors=False,
ref_dmdtda=ref_dmdtda, err_ref_dmdtda=0.000001,
maxiter_mu_star=2)
maxiter=2)
# test that fallback function works as expected if ignore_error=True and
# if the first guess can improve (but not enough)
model_fallback = run_dynamic_mu_star_calibration(
gdir, max_mu_star=1000.,
run_function=dynamic_mu_star_run_with_dynamic_spinup,
kwargs_run_function={'minimise_for': minimise_for,
'precision_percent_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'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_dyn_spinup': precision_percent,
'precision_absolute_dyn_spinup': precision_absolute,
'precision_percent': precision_percent,
'precision_absolute': precision_absolute,
'do_inversion': do_inversion},
output_filesuffix='_dyn_mu_calib_spinup_inversion_error',
ignore_errors=True,
ref_dmdtda=ref_dmdtda, err_ref_dmdtda=0.000001,
maxiter_mu_star=2)
maxiter=2)
assert isinstance(model_fallback, oggm.core.flowline.FluxBasedModel)
assert gdir.get_diagnostics()['used_spinup_option'] == \
'dynamic mu_star calibration (part success)'
Expand All @@ -3944,20 +3944,20 @@ def test_run_dynamic_mu_star_calibration_with_dynamic_spinup(self,
gdir, max_mu_star=1000.,
run_function=dynamic_mu_star_run_with_dynamic_spinup,
kwargs_run_function={'minimise_for': minimise_for,
'precision_percent_dyn_spinup': 0.1,
'precision_absolute_dyn_spinup': 0.0001,
'maxiter_dyn_spinup': 2,
'precision_percent': 0.1,
'precision_absolute': 0.0001,
'maxiter': 2,
'do_inversion': do_inversion},
fallback_function=dynamic_mu_star_run_with_dynamic_spinup_fallback,
kwargs_fallback_function={'minimise_for': minimise_for,
'precision_percent_dyn_spinup': 0.1,
'precision_absolute_dyn_spinup': 0.0001,
'maxiter_dyn_spinup': 2,
'precision_percent': 0.1,
'precision_absolute': 0.0001,
'maxiter': 2,
'do_inversion': do_inversion},
output_filesuffix='_dyn_mu_calib_spinup_inversion_error',
ignore_errors=True,
ref_dmdtda=ref_dmdtda, err_ref_dmdtda=0.000001,
maxiter_mu_star=2)
maxiter=2)
assert isinstance(model_fallback, oggm.core.flowline.FluxBasedModel)
assert gdir.get_diagnostics()['used_spinup_option'] == \
'fixed geometry spinup'
Expand Down Expand Up @@ -4151,7 +4151,7 @@ def test_run_dynamic_mu_star_calibration_without_dynamic_spinup(self):
output_filesuffix='_dyn_mu_calib_error',
ignore_errors=False,
ref_dmdtda=ref_dmdtda, err_ref_dmdtda=0.000001,
maxiter_mu_star=2)
maxiter=2)
# test that fallback function works as expected if ignore_error=True and
# if the first guess can improve (but not enough)
model_fallback = run_dynamic_mu_star_calibration(
Expand All @@ -4161,7 +4161,7 @@ def test_run_dynamic_mu_star_calibration_without_dynamic_spinup(self):
output_filesuffix='_dyn_mu_calib_spinup_inversion_error',
ignore_errors=True,
ref_dmdtda=ref_dmdtda, err_ref_dmdtda=0.000001,
maxiter_mu_star=2)
maxiter=2)
assert isinstance(model_fallback, oggm.core.flowline.FluxBasedModel)
assert gdir.get_diagnostics()['used_spinup_option'] == \
'dynamic mu_star calibration (part success)'
Expand All @@ -4176,7 +4176,7 @@ def test_run_dynamic_mu_star_calibration_without_dynamic_spinup(self):
output_filesuffix='_dyn_mu_calib_error',
ignore_errors=True,
ref_dmdtda=ref_dmdtda, err_ref_dmdtda=0.000001,
maxiter_mu_star=2)
maxiter=2)
assert isinstance(model_fallback, oggm.core.flowline.FluxBasedModel)
assert gdir.get_diagnostics()['used_spinup_option'] == 'no spinup'

Expand Down
2 changes: 1 addition & 1 deletion oggm/tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -1247,7 +1247,7 @@ def test_elev_bands_and_spinup_run(self):
# Around RGI date they are close
# have to exclude rgi_id 'RGI60-11.00719_d02', because no dmdtda data
assert_allclose(dss.sel(time=2004).area[1:], dse.sel(time=2004).area[1:], rtol=0.01)
assert_allclose(dss.sel(time=2004).length[1:], dse.sel(time=2004).length[1:], atol=805)
assert_allclose(dss.sel(time=2004).length[1:], dse.sel(time=2004).length[1:], atol=940)
assert_allclose(dss.sel(time=2004).volume[1:], dse.sel(time=2004).volume[1:], rtol=0.21)

# Over the period they are... close enough
Expand Down
1 change: 0 additions & 1 deletion oggm/tests/test_workflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,6 @@ def test_calibrate_inversion_from_consensus(self):
gdirs, ignore_missing=True,
volume_m3_reference=user_provided_volume_m3)

df = df.dropna()
np.testing.assert_allclose(user_provided_volume_m3,
df.vol_oggm_m3.sum(),
rtol=0.01)
Expand Down
7 changes: 6 additions & 1 deletion oggm/workflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -740,7 +740,12 @@ def compute_vol(x):
filter_inversion_output=filter_inversion_output)
odf = df.copy()
odf['oggm'] = execute_entity_task(tasks.get_inversion_volume, gdirs)
return odf.dropna()
# if the user provides a glacier volume all glaciers are considered,
# dropna() below exclude glaciers where no ITMIX volume is available
if volume_m3_reference is None:
return odf.dropna()
else:
return odf

def to_minimize(x):
log.workflow('Consensus estimate optimisation with '
Expand Down