diff --git a/docs/whats-new.rst b/docs/whats-new.rst index bff58256a..a2b55b57e 100644 --- a/docs/whats-new.rst +++ b/docs/whats-new.rst @@ -6,10 +6,18 @@ Version history v1.6.2 (unreleased) ------------------- +Enhancements +~~~~~~~~~~~~ + +- There is now a possibility for initializing a elevation-band flowline using + external thickness data and conduct a dynamic run with it (:pull:`1658`). + Bug fixes ~~~~~~~~~ -- The binned variables in the elevation band flowlines did not used the - glacier mask when preserving the total values (:pull:`1661`). + +- The binned variables in the elevation band flowlines did not use the + glacier mask when preserving the total values. This is a bad + bug that is now fixed (:pull:`1661`). By `Patrick Schmitt `_ diff --git a/oggm/core/flowline.py b/oggm/core/flowline.py index 3ec746e6b..98167b67d 100644 --- a/oggm/core/flowline.py +++ b/oggm/core/flowline.py @@ -2996,7 +2996,8 @@ def calving_glacier_downstream_line(line, n_points): @entity_task(log, writes=['model_flowlines']) -def init_present_time_glacier(gdir, filesuffix=''): +def init_present_time_glacier(gdir, filesuffix='', + use_binned_thickness_data=False): """Merges data from preprocessing tasks. First task after inversion! This updates the `mode_flowlines` file and creates a stand-alone numerical @@ -3007,9 +3008,14 @@ def init_present_time_glacier(gdir, filesuffix=''): gdir : :py:class:`oggm.GlacierDirectory` the glacier directory to process filesuffix : str - append a suffix to the model_flowlines filename (e.g. useful for - dynamic melt_f calibration including an inversion, so the original - model_flowlines are not changed). + append a suffix to the model_flowlines filename (e.g. useful for + dynamic melt_f calibration including an inversion, so the original + model_flowlines are not changed). + use_binned_thickness_data : bool or str + if you want to use thickness data, which was binned to the elevation + band flowlines with tasks.elevation_band_flowine and + tasks.fixed_dx_elevation_band_flowline, you can provide the name of the + data here to create a flowline for a dynamic model run """ # Some vars @@ -3026,20 +3032,53 @@ def init_present_time_glacier(gdir, filesuffix=''): # Get the data to make the model flowlines line = cl.line - section = inv['volume'] / (cl.dx * map_dx) surface_h = cl.surface_h - bed_h = surface_h - inv['thick'] widths_m = cl.widths * map_dx - assert np.all(widths_m > 0) - bed_shape = 4 * inv['thick'] / (cl.widths * map_dx) ** 2 - lambdas = inv['thick'] * np.NaN - lambdas[inv['is_trapezoid']] = def_lambda - lambdas[inv['is_rectangular']] = 0. + if not use_binned_thickness_data: + # classical initialisation after the inversion + section = inv['volume'] / (cl.dx * map_dx) + bed_h = surface_h - inv['thick'] + + bed_shape = 4 * inv['thick'] / widths_m ** 2 + + lambdas = inv['thick'] * np.NaN + lambdas[inv['is_trapezoid']] = def_lambda + lambdas[inv['is_rectangular']] = 0. - # Where the flux and the thickness is zero we just assume trapezoid: - lambdas[bed_shape == 0] = def_lambda + # Where the flux and the thickness is zero we just assume trapezoid: + lambdas[bed_shape == 0] = def_lambda + + else: + # here we use binned thickness data for the initialisation + elev_fl = pd.read_csv( + gdir.get_filepath('elevation_band_flowline', + filesuffix='_fixed_dx'), index_col=0) + assert np.allclose(widths_m, elev_fl['widths_m']) + elev_fl_thick = elev_fl[use_binned_thickness_data].values + section = elev_fl_thick * widths_m + lambdas = np.ones(len(section)) * def_lambda + + # for trapezoidal the calculation of thickness results in quadratic + # equation, we only keep solution which results in a positive w0 + # value (-> w/lambda >= h), + # it still could happen that we would need a negative w0 if the + # section is too large, in those cases we get a negative value + # inside sqrt (we ignore the RuntimeWarning) -> if this happens we + # use rectangular shape with original thickness + with np.errstate(invalid='ignore'): + thick = ((2 * widths_m - + np.sqrt(4 * widths_m ** 2 - + 4 * lambdas * 2 * section)) / + (2 * lambdas)) + nan_thick = np.isnan(thick) + thick[nan_thick] = elev_fl_thick[nan_thick] + lambdas[nan_thick] = 0 + + # finally the glacier bed and other stuff + bed_h = surface_h - thick + bed_shape = 4 * thick / widths_m ** 2 if not gdir.is_tidewater and inv['is_last']: # for valley glaciers, simply add the downstream line, depending on diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index 2a6ed1648..34966d0c6 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -17,7 +17,7 @@ import xarray as xr from oggm import utils, workflow, tasks, cfg from oggm.core import climate, inversion, centerlines -from oggm.shop import gcm_climate +from oggm.shop import gcm_climate, bedtopo from oggm.cfg import SEC_IN_YEAR, SEC_IN_MONTH from oggm.utils import get_demo_file from oggm.exceptions import InvalidParamsError, InvalidWorkflowError @@ -153,6 +153,67 @@ def test_init_present_time_glacier(self, hef_gdir, downstream_line_shape): with pytest.raises(InvalidParamsError): init_present_time_glacier(gdir) + def test_init_present_time_glacier_obs_thick(self, hef_elev_gdir, + monkeypatch): + + gdir = hef_elev_gdir + + # need to change rgi_id, which is needed to be different in other tests + # when comparing to centerlines + gdir.rgi_id = 'RGI60-11.00897' + + # add some thickness data + ft = utils.get_demo_file('RGI60-11.00897_thickness.tif') + monkeypatch.setattr(utils, 'file_downloader', lambda x: ft) + bedtopo.add_consensus_thickness(gdir) + vn = 'consensus_ice_thickness' + centerlines.elevation_band_flowline(gdir, bin_variables=[vn]) + centerlines.fixed_dx_elevation_band_flowline(gdir, + bin_variables=[vn]) + + tasks.init_present_time_glacier(gdir, filesuffix='_consensus', + use_binned_thickness_data=vn) + fl_consensus = gdir.read_pickle('model_flowlines', + filesuffix='_consensus')[0] + + # check that resulting flowline has the same volume as observation + cdf = pd.read_hdf(utils.get_demo_file('rgi62_itmix_df.h5')) + ref_vol = cdf.loc[gdir.rgi_id].vol_itmix_m3 + np.testing.assert_allclose(fl_consensus.volume_m3, ref_vol) + + # should be trapezoid where ice + assert np.all(fl_consensus.is_trapezoid[fl_consensus.thick > 0]) + + # test that we can use fl in an dynamic model run without an error + mb = LinearMassBalance(3000.) + model_ref = FluxBasedModel(gdir.read_pickle('model_flowlines'), + mb_model=mb) + model_ref.run_until(100) + model_consensus = FluxBasedModel([fl_consensus], mb_model=mb) + model_consensus.run_until(100) + np.testing.assert_allclose(model_ref.volume_km3, + model_consensus.volume_km3, + atol=0.02) + + # test that if w0<0 it is converted to rectangular + # set some thickness to very large values to force it + df_fixed_dx = pd.read_csv(gdir.get_filepath('elevation_band_flowline', + filesuffix='_fixed_dx')) + new_thick = df_fixed_dx['consensus_ice_thickness'] + new_thick[-10:] = new_thick[-10:] + 1000 + df_fixed_dx['consensus_ice_thickness'] = new_thick + ref_vol_rect = np.sum(df_fixed_dx['area_m2'] * new_thick) + df_fixed_dx.to_csv(gdir.get_filepath('elevation_band_flowline', + filesuffix='_fixed_dx')) + + tasks.init_present_time_glacier(gdir, filesuffix='_consensus_rect', + use_binned_thickness_data=vn) + fl_consensus_rect = gdir.read_pickle('model_flowlines', + filesuffix='_consensus_rect')[0] + + np.testing.assert_allclose(fl_consensus_rect.volume_m3, ref_vol_rect) + assert np.sum(fl_consensus_rect.is_rectangular) == 10 + def test_present_time_glacier_massbalance(self, hef_gdir): gdir = hef_gdir