Compute smoother centerlines

This notebook describes how to compute centerlines with OGGM and write them to disk. It is meant for users who are mostly interested in the centerlines, not so much the rest of the OGGM model.

We use an example of a user-provided glacier inventory and DEM (thanks to Liss Andreassen for providing the data).

import rioxarray as rioxr
import geopandas as gpd
import matplotlib.pyplot as plt

from oggm import cfg, utils, workflow, tasks
ModuleNotFoundError                       Traceback (most recent call last)
Input In [1], in <module>
----> 1 import rioxarray as rioxr
      2 import geopandas as gpd
      3 import matplotlib.pyplot as plt

ModuleNotFoundError: No module named 'rioxarray'

Data preparation

Download the demo data. This is a subset of a regional glacier inventory and DEM in Norway:

fpath_inventory = utils.file_downloader('')
fpath_dem = utils.file_downloader('')

Read the data and plot it:

inventory = gpd.read_file(fpath_inventory)
dem = rioxr.open_rasterio(fpath_dem)

f, ax = plt.subplots(figsize=(9, 9))
dem.plot(ax=ax, cmap='terrain', vmin=0);
inventory.plot(ax=ax, edgecolor='k', facecolor='C1');

The resolution of the DEM is 10m:

print((dem.x[1] - dem.x[0]).item())

In this inventory, one geometry has a topological error (the figure of eight where the outlines touch):


Let’s correct it:

inventory.loc[~inventory.is_valid, 'geometry'] = inventory.loc[~inventory.is_valid].buffer(0)

A final preparation step is to convert the format of the inventory to a file which resembles the RGI (see use_your_own_inventory.ipynb):

# We keep the original ID for later reference
gdf = utils.cook_rgidf(inventory, o1_region='08', assign_column_values={'breID':'breID'})

Compute the centerlines

We use the standard OGGM procedure for this:


# Parameters
cfg.PARAMS['use_multiprocessing'] = True  # this is often a good idea
cfg.PARAMS['use_rgi_area'] = False  # this is required for user-defined inventories
cfg.PARAMS['use_intersects'] = False  # we don't care about intersects for centerlines
cfg.PARAMS['border'] = 10  # no need to make a large map

# Optional: change the grid resolution
# E.g. fixed grid spacing
# cfg.PARAMS['grid_dx_method'] = 'fixed'
# cfg.PARAMS['fixed_dx'] = 10
# Or variable but twice higher than default 
cfg.PARAMS['grid_dx_method'] = 'square'
cfg.PARAMS['d1'] = 7  # (default is 14)
cfg.PARAMS['d2'] = 5  # (default is 10)
cfg.PARAMS['dmax'] = 200  # (default is 100)

# Tell OGGM to use our user DEM (important!)
cfg.PATHS['dem_file'] = fpath_dem

# Where to work
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='NORWAY_CENTERLINES', reset=True)

Now the workflow:

gdirs = workflow.init_glacier_directories(gdf)

workflow.execute_entity_task(tasks.define_glacier_region, gdirs, source='USER');  # Use the user DEM

workflow.execute_entity_task(tasks.glacier_masks, gdirs);
workflow.execute_entity_task(tasks.compute_centerlines, gdirs);

Note: the default in OGGM is to use a grid size of varying resolution for each glacier. I think it makes sense in many cases, but you may prefer to use the native resolution of your DEM. You can do so by commenting / un-commenting the options above.

Write the data to a shapefile with optional smoothing

The relevant task is “write_centerlines_to_shape”, which writes everything to a shapefile:

from oggm.utils import write_centerlines_to_shape

write_centerlines_to_shape(gdirs,  # The glaciers to process
                           path='outputs/Norway_Centerlines.shp',  # The output file
                           to_tar=False,  # set to True to put everything into one single tar file
                 ,  # Write into the projection of the original inventory
                           keep_main_only=True,  # Write only the main flowline and discard the tributaries

Let’s have a look at the output:

cls_default = gpd.read_file('outputs/Norway_Centerlines.shp')
cls_default['breID'] = gdf['breID']  # This only works this way because we have one centerline per glacier!


LE_SEGMENT is the length of the centerline in meters. The RGI “IDs” are fake (OGGM needs them) but the breID are real. Lets use them as index for the file:

cls_default = cls_default.set_index('breID')
orig_inventory = inventory.set_index('breID')

Now we can plot an example:

sel_breID = 1189  # 5570

f, ax = plt.subplots(figsize=(9, 4))
orig_inventory.loc[[sel_breID]].plot(ax=ax, facecolor='lightblue');

What can we see?

-the centerline does not end exactly at the glacier outline -the line seems “crooked”, it has sudden turns

Both effects are due to the algorithm we use to compute the centerlines (Kienholz et al., (2014)), which works on the underlying glacier grid. Each vertice (point) in the line corresponds to the center of the grid point.

We have implemented a few new options in OGGM v1.6, which allow to circumvent these limitations. We illustrate them here:

write_centerlines_to_shape(gdirs,  # The glaciers to process
                           path='outputs/Norway_Centerlines_smooth.shp',  # The output file
                           to_tar=False,  # set to True to put everything into one single tar file
                 ,  # Write into the projection of the original inventory
                           keep_main_only=True,  # Write only the main flowline and discard the tributaries
                           ensure_exterior_match=True,  # NEW! Ensure that the lines are touching the outlines
                           simplify_line=0.5,  # NEW! this option reduces the number of vertices along the line
                           corner_cutting=5,  # BEW! this then augments the number of vertices again

The simplify_line and corner_cutting options are cosmetic and subjective. The former will simplify the line, by making it look less edgy but also less precise, while the latter then “smoothes” it. Users may try different combinations to see their effect (see the documentation).

cls_smooth = gpd.read_file('outputs/Norway_Centerlines_smooth.shp')
cls_smooth['breID'] = gdf['breID']
cls_smooth = cls_smooth.set_index('breID')
sel_breID = 1189

f, ax = plt.subplots(figsize=(9, 4))
orig_inventory.loc[[sel_breID]].plot(ax=ax, facecolor='lightblue');
cls_smooth.loc[[sel_breID]].plot(ax=ax, color='C3');

Final remarks

While the centerline algorithm is quite robust, the results will vary as a function of the resolution of the underlying grid, and the smoothing options. After trying a little, it seems difficult to find a setting which works “best” in all circumstances, and we encourage users to try several options and see what they prefer. The option likely to have the most impact (assuming smoothing with (0.5, 5) is the underlying grid resolution.

What’s next?