OGGM flowlines: where are they?

In this notebook we show how to access the OGGM flowlines location before, during, and after a run.

Some of the code shown here will make it to the OGGM codebase eventually.

from oggm import cfg, utils, workflow, tasks, graphics
from oggm.core import flowline
import salem
import xarray as xr
import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
cfg.initialize(logging_level='WARNING')
2021-09-15 10:43:21: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2021-09-15 10:43:21: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2021-09-15 10:43:21: oggm.cfg: Multiprocessing: using all available processors (N=2)

Get ready

# Where to store the data 
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-flowlines', reset=True)
# Which glaciers?
rgi_ids = ['RGI60-11.00897']
# We start from prepro level 3 with all data ready
gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_border=40)
gdir = gdirs[0]
gdir
2021-09-15 10:43:22: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2021-09-15 10:43:22: oggm.workflow: Execute entity task gdir_from_prepro on 1 glaciers
<oggm.GlacierDirectory>
  RGI id: RGI60-11.00897
  Region: 11: Central Europe
  Subregion: 11-01: Alps                            
  Name: Hintereisferner
  Glacier type: Glacier
  Terminus type: Land-terminating
  Status: Glacier or ice cap
  Area: 8.036 km2
  Lon, Lat: (10.7584, 46.8003)
  Grid (nx, ny): (199, 154)
  Grid (dx, dy): (50.0, -50.0)

Where is the terminus of the RGI glacier?

There are several ways to get the terminus, depending on what you want. They are also not necessarily exact same:

Terminus as the lowest point on the glacier

# Get the topo data and the glacier mask
with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    topo = ds.topo
    # Glacier outline raster
    mask = ds.glacier_ext
topo.plot();
../_images/where_are_the_flowlines_10_0.png
topo_ext = topo.where(mask==1)
topo_ext.plot();
../_images/where_are_the_flowlines_11_0.png
# Get the terminus
terminus = topo_ext.where(topo_ext==topo_ext.min(), drop=True)
# Project its coordinates from the local UTM to WGS-84
t_lon, t_lat = salem.transform_proj(gdir.grid.proj, 'EPSG:4326', terminus.x[0], terminus.y[0])
print('lon, lat:', t_lon, t_lat)
print('google link:', f'https://www.google.com/maps/place/{t_lat},{t_lon}')
lon, lat: 10.802746863388563 46.818900914354266
google link: https://www.google.com/maps/place/46.818900914354266,10.802746863388563

Terminus as the lowest point on the main centerline

# Get the centerlines
cls = gdir.read_pickle('centerlines')
# Get the coord of the last point of the main centerline
cl = cls[-1]
i, j = cl.line.coords[-1]
# These coords are in glacier grid coordinates. Let's convert them to lon, lat:
t_lon, t_lat = gdir.grid.ij_to_crs(i, j, crs='EPSG:4326')
print('lon, lat:', t_lon, t_lat)
print('google link:', f'https://www.google.com/maps/place/{t_lat},{t_lon}')
lon, lat: 10.802746863955196 46.818901155282155
google link: https://www.google.com/maps/place/46.818901155282155,10.802746863955196

Terminus as the lowest point on the main flowline

“centerline” in the OGGM jargon is not the same as “flowline”. Flowlines have a fixed dx and their terminus is not necessarily exact on the glacier outline. Code-wise it’s very similar though:

# Get the flowlines
cls = gdir.read_pickle('inversion_flowlines')
# Get the coord of the last point of the main centerline
cl = cls[-1]
i, j = cl.line.coords[-1]
# These coords are in glacier grid coordinates. Let's convert them to lon, lat:
t_lon, t_lat = gdir.grid.ij_to_crs(i, j, crs='EPSG:4326')
print('lon, lat:', t_lon, t_lat)
print('google link:', f'https://www.google.com/maps/place/{t_lat},{t_lon}')
lon, lat: 10.802746777546085 46.818796056851475
google link: https://www.google.com/maps/place/46.818796056851475,10.802746777546085

Bonus: convert the centerlines to a shapefile

utils.write_centerlines_to_shape(gdirs, path='centerlines.shp')
2021-09-15 10:43:24: oggm.utils: Applying global task write_centerlines_to_shape on 1 glaciers
2021-09-15 10:43:24: oggm.utils: write_centerlines_to_shape on centerlines.shp ...
2021-09-15 10:43:24: oggm.workflow: Execute entity task _get_centerline_lonlat on 1 glaciers
/usr/local/pyenv/versions/3.8.12/lib/python3.8/site-packages/pyproj/crs/crs.py:131: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  in_crs_string = _prepare_from_proj_string(in_crs_string)
sh = gpd.read_file('centerlines.shp')
sh.plot();
../_images/where_are_the_flowlines_20_0.png

Remember: the “centerlines” are not the same things as “flowlines” in OGGM. The later objects undergo further quality checks, such as the impossibility for ice to “climb”, i.e. have negative slopes. The flowlines are therefore sometimes shorter than the centerlines:

utils.write_centerlines_to_shape(gdirs, path='flowlines.shp', flowlines_output=True)
sh = gpd.read_file('flowlines.shp')
sh.plot();
2021-09-15 10:43:24: oggm.utils: Applying global task write_centerlines_to_shape on 1 glaciers
2021-09-15 10:43:24: oggm.utils: write_centerlines_to_shape on flowlines.shp ...
2021-09-15 10:43:24: oggm.workflow: Execute entity task _get_centerline_lonlat on 1 glaciers
/usr/local/pyenv/versions/3.8.12/lib/python3.8/site-packages/pyproj/crs/crs.py:131: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  in_crs_string = _prepare_from_proj_string(in_crs_string)
../_images/where_are_the_flowlines_22_4.png

Flowline geometry after a run with FileModel

Let’s do a run first:

tasks.init_present_time_glacier(gdir)
tasks.run_constant_climate(gdir, nyears=100, y0=2000);

We use a FileModel to read the model output:

fmod = flowline.FileModel(gdir.get_filepath('model_geometry'))

A FileModel behaves like a OGGM’s FlowlineModel:

fmod.run_until(0)  # Point the file model to year 0 in the output
graphics.plot_modeloutput_map(gdir, model=fmod)  # plot it
../_images/where_are_the_flowlines_29_0.png
fmod.run_until(100)  # Point the file model to year 100 in the output
graphics.plot_modeloutput_map(gdir, model=fmod)  # plot it
../_images/where_are_the_flowlines_30_0.png
# Bonus - get back to e.g. the volume timeseries
fmod.volume_km3_ts().plot();
../_images/where_are_the_flowlines_31_0.png

OK, now create a table of the main flowline’s grid points location and bed altitude (this does not change with time):

fl = fmod.fls[-1]  # Main flowline
i, j = fl.line.xy  # xy flowline on grid
lons, lats = gdir.grid.ij_to_crs(i, j, crs='EPSG:4326')  # to WGS84

df_coords = pd.DataFrame(index=fl.dis_on_line*gdir.grid.dx)
df_coords.index.name = 'Distance along flowline'
df_coords['lon'] = lons
df_coords['lat'] = lats
df_coords['bed_elevation'] = fl.bed_h
df_coords.plot(x='lon', y='lat');
../_images/where_are_the_flowlines_33_0.png
df_coords['bed_elevation'].plot();
../_images/where_are_the_flowlines_34_0.png

Now store a time varying array of ice thickness, surface elevation along this line:

years = np.arange(0, 101)
df_thick = pd.DataFrame(index=df_coords.index, columns=years, dtype=np.float64)
df_surf_h = pd.DataFrame(index=df_coords.index, columns=years, dtype=np.float64)
df_bed_h = pd.DataFrame()
for year in years:
    fmod.run_until(year)
    fl = fmod.fls[-1]
    df_thick[year] = fl.thick
    df_surf_h[year] = fl.surface_h
df_thick[[0, 50, 100]].plot();
plt.title('Ice thickness at three points in time');
../_images/where_are_the_flowlines_37_0.png
f, ax = plt.subplots()
df_surf_h[[0, 50, 100]].plot(ax=ax);
df_coords['bed_elevation'].plot(ax=ax, color='k');
plt.title('Glacier elevation at three points in time');
../_images/where_are_the_flowlines_38_0.png

Location of the terminus over time

Let’s find the indices where the terminus is (i.e. the last point where ice is thicker than 1m), and link these to the lon, lat positions along the flowlines.

The first method uses fancy pandas functions but may be more cryptic for less experienced pandas users:

# Nice trick from https://stackoverflow.com/questions/34384349/find-index-of-last-true-value-in-pandas-series-or-dataframe
dis_term = (df_thick > 1)[::-1].idxmax()

# Select the terminus coordinates at these locations
loc_over_time = df_coords.loc[dis_term].set_index(dis_term.index)

# Plot them over time
loc_over_time.plot.scatter(x='lon', y='lat', c=loc_over_time.index, colormap='viridis');
plt.title('Location of the terminus over time');
../_images/where_are_the_flowlines_41_0.png
# Plot them on a google image - you need an API key for this
# api_key = ''
# from motionless import DecoratedMap, LatLonMarker
# dmap = DecoratedMap(maptype='satellite', key=api_key)
# for y in [0, 20, 40, 60, 80, 100]:
#     tmp = loc_over_time.loc[y]
#     dmap.add_marker(LatLonMarker(tmp.lat, tmp.lon, ))
# print(dmap.generate_url())

And now, method 2: less fancy but maybe easier to read?

for yr in [0, 20, 40, 60, 80, 100]:
    # Find the last index of the terminus
    p_term = np.nonzero(df_thick[yr].values > 1)[0][-1]
    # Print the location of the terminus
    print(f'Terminus pos at year {yr}', df_coords.iloc[p_term][['lon', 'lat']].values)
Terminus pos at year 0 [10.80274678 46.81879606]
Terminus pos at year 20 [10.7936729  46.81537664]
Terminus pos at year 40 [10.77756361 46.80792269]
Terminus pos at year 60 [10.76641209 46.79531908]
Terminus pos at year 80 [10.75236938 46.79236233]
Terminus pos at year 100 [10.75236938 46.79236233]

Comments on “elevation band flowlines”

If you use elevation band flowlines, the location of the flowlines is not known: indeed, the glacier is an even more simplified representation of the real world one. In this case, if you are interested in tracking the terminus position, you may need to use tricks, such as using the retreat from the terminus with time, or similar.

What’s next?