Computing the centerlines is the first task to run after the initialisation of the local glacier directories and of the local topography.
Our algorithm is an implementation of the procedure described by Kienholz et al., (2014). Appart from some minor changes (mostly the choice of certain parameters), we stayed close to the original algorithm.
The relevant task is
In : graphics.plot_centerlines(gdir)
The glacier has a major centerline (the longest one), and
tributaries (in this case two ). The
Centerline objects are stored
as a list, the last one being
the major one. Navigation between inflows (can be more than one) and
outflow (only one or none) is facilitated by the
In : fls = gdir.read_pickle('centerlines') In : fls # a Centerline object Out: <oggm.core.centerlines.Centerline at 0x7f09dc5546d8> # make sure the first flowline realy flows into the major one: In : assert fls.flows_to is fls[-1]
At this stage, the centerline coordinates are still defined on the original
grid, and they are not considered as “flowlines” by OGGM. A rather simple task
tasks.initialize_flowlines()) converts them to flowlines which
now have a regular coordinate spacing along the flowline (which they will
keep for the rest of the workflow). The tail of the tributaries are cut
according to a distance threshold rule:
In : graphics.plot_centerlines(gdir, use_flowlines=True)
For the glacier to be able to grow we need to determine the flowlines
downstream of the current glacier geometry. This is done by the
In : graphics.plot_centerlines(gdir, use_flowlines=True, add_downstream=True)
The downsteam lines area also computed using a routing algorithm minimizing the distance to cover and upward slopes.
Each flowline has it’s own “catchment area”. These areas are computed using similar flow routing methods as the one used for determining the flowlines. Their role is to attribute each glacier pixel to the right tributory (this will also influence the later computation of the glacier widths).
In : tasks.catchment_area(gdir) In : graphics.plot_catchment_areas(gdir)
Finally, the glacier widths are computed in two steps.
First, we compute the geometrical width at each grid point. The width is drawn from the intersection of a line normal to the flowline and either the glacier or the catchment outlines (when there are tributaries):
In : tasks.catchment_width_geom(gdir) I am densified (external_values, 10 elements) In : graphics.plot_catchment_width(gdir)
Then, these geometrical widths are corrected so that the altitude-area
distribution of the “flowline-glacier” is as close as possible as the actual
distribution of the glacier using its full 2D geometry. This job is done by
In : tasks.catchment_width_correction(gdir) In : graphics.plot_catchment_width(gdir, corrected=True)
Note that a perfect distribution is not possible since the sample size is not the same between the “1.5D” and the 2D representation of the glacier. OGGM deals with this by iteratively search for an altidute bin size which ensures that both representations have at least one element for each bin.
Shared setup for these examples:
import os import geopandas as gpd import oggm from oggm import cfg, tasks from oggm.utils import get_demo_file cfg.initialize() cfg.set_intersects_db(get_demo_file('rgi_intersect_oetztal.shp')) cfg.PATHS['dem_file'] = get_demo_file('hef_srtm.tif') base_dir = os.path.join(oggm.gettempdir(), 'OGGM_docs', 'Flowlines') entity = gpd.read_file(get_demo_file('HEF_MajDivide.shp')).iloc gdir = oggm.GlacierDirectory(entity, base_dir=base_dir) tasks.define_glacier_region(gdir, entity=entity) tasks.glacier_masks(gdir) tasks.compute_centerlines(gdir) tasks.initialize_flowlines(gdir) tasks.compute_downstream_line(gdir)