OGGM flowlines: where are they?

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')

2024-08-25 21:30:17: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:30:17: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:30:17: oggm.cfg: Multiprocessing: using all available processors (N=4)

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
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/centerlines/W5E5/'
gdirs = workflow.init_glacier_directories(rgi_ids, from_prepro_level=3, prepro_base_url=base_url, prepro_border=80)
gdir = gdirs[0]
gdir

2024-08-25 21:30:18: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2024-08-25 21:30:18: oggm.workflow: Execute entity tasks [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): (279, 234)
  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/e66f52c23f57d362a68ea6f1fec8aa15bb0e713745299ecc6216875f94fb2531.png

topo_ext = topo.where(mask==1)
topo_ext.plot();

../../_images/c50d1befa4ecac5112b69cff4730788da80ce18d36ae8740c28c50af0049006e.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#

output_dir = utils.mkdir('outputs')
utils.write_centerlines_to_shape(gdirs, path=f'{output_dir}/centerlines.shp')

2024-08-25 21:30:20: oggm.utils: Applying global task write_centerlines_to_shape on 1 glaciers
2024-08-25 21:30:20: oggm.utils: write_centerlines_to_shape on outputs/centerlines.shp ...
2024-08-25 21:30:20: oggm.workflow: Execute entity tasks [get_centerline_lonlat] on 1 glaciers

sh = gpd.read_file(f'{output_dir}/centerlines.shp')
sh.plot();

../../_images/f7c507ed87e818b065bf3ed7ecb37d792b4a6481e3f6b9281edbe111a5e6bdb7.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=f'{output_dir}/flowlines.shp', flowlines_output=True)
sh = gpd.read_file(f'{output_dir}/flowlines.shp')
sh.plot();

2024-08-25 21:30:20: oggm.utils: Applying global task write_centerlines_to_shape on 1 glaciers
2024-08-25 21:30:20: oggm.utils: write_centerlines_to_shape on outputs/flowlines.shp ...
2024-08-25 21:30:20: oggm.workflow: Execute entity tasks [get_centerline_lonlat] on 1 glaciers

../../_images/92b1702ff0d9f67a991213a1ba2b875afe3a0a5812f9b04f9c5c050f75b76b3b.png

Flowline geometry after a run: with `FileModel`#

This method uses another way to output the model geometry (model_geometry files), which are less demanding in disk space but require more computations to restore the full geometry after a run (OGGM will help you to do that). This method is to be preferred to the above if you want to save disk space or if you don’t have yet access to flowline diagnostics files.

Let’s do a run first:

cfg.PARAMS['store_model_geometry'] = True  # We want to get back to it later

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

2024-08-25 21:30:23: oggm.cfg: PARAMS['store_model_geometry'] changed from `False` to `True`.

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/f4f0d73b101c94ed929c206ecb85a8ca144e782a4a174e392cc8fdcb0a111dd3.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/e6c96f4b838a7e4ec478cfd17e2ad08545c8e4935a62d08ca4e71b5846d8d6e6.png

# Bonus - get back to e.g. the volume timeseries
fmod.volume_km3_ts().plot();

../../_images/72844be3d7dcd6de7589eb8e81cff7e4cc5e643f181dde007d3df12a25f5ae9e.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', legend=False)
plt.ylabel('lat');

../../_images/839fe11123ef2a6f623e54d316e379159c2ddcba17714beb549368e657904a4e.png

df_coords['bed_elevation'].plot();

../../_images/a27c6afdb340ec154971369b7c67a54ab0c90aeb6f06cff0c2591490bee0596b.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')
plt.ylabel('Glacier thickness (m)')
plt.legend(title='Year');

../../_images/3c4adff2c83298af0973c9b2ce47f073d99963c3bf6bc7d6667eb00b766aa630.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')
plt.legend(title='Year');

../../_images/2f1588a26444020e7685ee9b4b396aef261a2f0824d07245629790dbea4ca634.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)
loc_over_time['Year'] = loc_over_time.index.astype(int)

# Plot them over time
loc_over_time.plot.scatter(x='lon', y='lat', c='Year', colormap='viridis');
plt.title('Location of the terminus over time');

../../_images/8b8af32adfde29a5576d714c49b5aa9bf403f394dd9198055a203feda7cc6e86.png

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.80055249 46.81784614]
Terminus pos at year 40 [10.79152134 46.8143506 ]
Terminus pos at year 60 [10.77663662 46.80728651]
Terminus pos at year 80 [10.76826544 46.7965916 ]
Terminus pos at year 100 [10.76154133 46.79236244]

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?#

return to the OGGM documentation
back to the table of contents

OGGM flowlines: where are they?

Contents

OGGM flowlines: where are they?#

Get ready#

Where is the terminus of the RGI glacier?#

Terminus as the lowest point on the glacier#

Terminus as the lowest point on the main centerline#

Terminus as the lowest point on the main flowline#

Bonus: convert the centerlines to a shapefile#

Flowline geometry after a run: with the new flowline diagnostics (new in v1.6.0!!)#

Location of the terminus over time#

Flowline geometry after a run: with `FileModel`#

Location of the terminus over time#

Comments on “elevation band flowlines”#

What’s next?#

OGGM flowlines: where are they?

Contents

OGGM flowlines: where are they?#

Get ready#

Where is the terminus of the RGI glacier?#

Terminus as the lowest point on the glacier#

Terminus as the lowest point on the main centerline#

Terminus as the lowest point on the main flowline#

Bonus: convert the centerlines to a shapefile#

Flowline geometry after a run: with the new flowline diagnostics (new in v1.6.0!!)#

Location of the terminus over time#

Flowline geometry after a run: with FileModel#

Location of the terminus over time#

Comments on “elevation band flowlines”#

What’s next?#

Flowline geometry after a run: with `FileModel`#