Skip to content

Commit

Permalink
Initialise dynamic flowline with observed thickness (#1658)
Browse files Browse the repository at this point in the history 
* add possibility for initialisation of a dynamic flowline with observed thickness data

* added whats-new

---------

Co-authored-by: Fabien Maussion <[email protected]>
  • Loading branch information
@pat-schmitt @fmaussion
pat-schmitt and fmaussion authored Nov 9, 2023
1 parent 365eef8 commit d40038b
Show file tree
Hide file tree
Showing 3 changed files with 124 additions and 16 deletions.
12 changes: 10 additions & 2 deletions docs/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <https://github.com/pat-schmitt>`_


Expand Down
65 changes: 52 additions & 13 deletions oggm/core/flowline.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
63 changes: 62 additions & 1 deletion oggm/tests/test_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit d40038b

Please sign in to comment.