RGI-TOPO for RGI 7.0#

OGGM was used to generate the topography data used to compute the topographical attributes and the centerlines products for RGI v7.0.

Here we show how to access this data from OGGM.

Input parameters#

This notebook can be run as a script with parameters using papermill, but it is not necessary. The following cell contains the parameters you can choose from:

# The RGI-id of the glaciers you want to look for
# Use the original shapefiles or the GLIMS viewer to check for the ID: https://www.glims.org/maps/glims
rgi_id = 'RGI2000-v7.0-G-01-06486'  # Denali

# The default is to test for all sources available for this glacier
# Set to a list of source names to override this
sources = None
# Where to write the plots. Default is in the current working directory
plot_dir = f'outputs/{rgi_id}'
# The RGI version to use
# V62 is an unofficial modification of V6 with only minor, backwards compatible modifications
prepro_rgi_version = 62
# Size of the map around the glacier. Currently only 10 and 40 are available
prepro_border = 10
# Degree of processing level.  Currently only 1 is available.
from_prepro_level = 1

Check input and set up#

# The sources can be given as parameters
if sources is not None and isinstance(sources, str):
    sources = sources.split(',')
# Plotting directory as well
if not plot_dir:
    plot_dir = './' + rgi_id
import os
plot_dir = os.path.abspath(plot_dir)
from oggm import cfg, utils, workflow, tasks, graphics, GlacierDirectory
import pandas as pd
import numpy as np
import xarray as xr
import rioxarray as rioxr
import geopandas as gpd
import salem
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
import itertools

from oggm.utils import DEM_SOURCES
from oggm.workflow import init_glacier_directories
# Make sure the plot directory exists
utils.mkdir(plot_dir);
# Use OGGM to download the data
cfg.initialize()
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-RGITOPO-RGI7', reset=True)
cfg.PARAMS['use_intersects'] = False
2026-08-07 07:53:19: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2026-08-07 07:53:19: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2026-08-07 07:53:19: oggm.cfg: Multiprocessing: using all available processors (N=4)
2026-08-07 07:53:19: oggm.cfg: PARAMS['use_intersects'] changed from `True` to `False`.

Download the data using OGGM utility functions#

Note that you could reach the same goal by downloading the data manually from

# URL of the preprocessed GDirs
gdir_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2026.1/all_dems/'
# We use OGGM to download the data
gdir = init_glacier_directories([rgi_id], from_prepro_level=1, prepro_border=10,  prepro_rgi_version='70G', prepro_base_url=gdir_url)[0]
2026-08-07 07:53:20: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2026-08-07 07:53:20: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
gdir
<oggm.GlacierDirectory>
  RGI id: RGI2000-v7.0-G-01-06486
  Region: 01: Alaska
  Subregion: 01-02: Alaska Range (Wrangell/Kilbuck)
  Glacier type: Glacier
  Terminus type: Not assigned
  Status: Glacier
  Area: 0.961122833004822 km2
  Lon, Lat: (-151.0094740399913, 63.061062)
  Grid (nx, ny): (70, 88)
  Grid (dx, dy): (24.0, -24.0)

Read the DEMs and store them all in a dataset#

if sources is None:
    sources = [src for src in os.listdir(gdir.dir) if src in utils.DEM_SOURCES]
print('RGI ID:', rgi_id)
print('Available DEM sources:', sources)
print('Plotting directory:', plot_dir)
RGI ID: RGI2000-v7.0-G-01-06486
Available DEM sources: ['ARCTICDEM', 'COPDEM90', 'ALASKA', 'DEM3', 'MAPZEN', 'TANDEM', 'COPDEM30', 'AW3D30', 'ASTER']
Plotting directory: /__w/tutorials/tutorials/notebooks/tutorials/outputs/RGI2000-v7.0-G-01-06486
# We use xarray to store the data
ods = xr.Dataset()
for src in sources:
    demfile = os.path.join(gdir.dir, src) + '/dem.tif'
    with rioxr.open_rasterio(demfile) as ds:
        data = ds.sel(band=1).load() * 1.
        ods[src] = data.where(data > -100, np.nan)

    sy, sx = np.gradient(ods[src], gdir.grid.dx, gdir.grid.dx)
    ods[src + '_slope'] = ('y', 'x'),  np.arctan(np.sqrt(sy**2 + sx**2))

with rioxr.open_rasterio(gdir.get_filepath('glacier_mask')) as ds:
    ods['mask'] = ds.sel(band=1).load()
# Decide on the number of plots and figure size
ns = len(sources)
x_size = 12
n_cols = 3
n_rows = -(-ns // n_cols)
y_size = x_size / n_cols * n_rows

Raw topography data#

smap = salem.graphics.Map(gdir.grid, countries=False)
smap.set_shapefile(gdir.read_shapefile('outlines'))
smap.set_plot_params(cmap='topo')
smap.set_lonlat_contours(add_tick_labels=False)
smap.set_plot_params(vmin=np.nanquantile([ods[s].min() for s in sources], 0.25),
                     vmax=np.nanquantile([ods[s].max() for s in sources], 0.75))

fig = plt.figure(figsize=(x_size, y_size))
grid = AxesGrid(fig, 111,
                nrows_ncols=(n_rows, n_cols),
                axes_pad=0.7,
                cbar_mode='each',
                cbar_location='right',
                cbar_pad=0.1
                )

for i, s in enumerate(sources):
    data = ods[s]
    smap.set_data(data)
    ax = grid[i]
    smap.visualize(ax=ax, addcbar=False, title=s)
    if np.isnan(data).all():
        grid[i].cax.remove()
        continue
    cax = grid.cbar_axes[i]
    smap.colorbarbase(cax)

# take care of uneven grids
if ax != grid[-1] and not grid[-1].title.get_text():
    grid[-1].remove()
    grid[-1].cax.remove()
if ax != grid[-2] and not grid[-2].title.get_text():
    grid[-2].remove()
    grid[-2].cax.remove()

plt.savefig(os.path.join(plot_dir, 'dem_topo_color.png'), dpi=150, bbox_inches='tight')
../../_images/d7d5af95e2c6477d01d01920ecae4d95ff3704533ff032af702149400e317e3a.png

Shaded relief#

fig = plt.figure(figsize=(x_size, y_size))
grid = AxesGrid(fig, 111,
                nrows_ncols=(n_rows, n_cols),
                axes_pad=0.7,
                cbar_location='right',
                cbar_pad=0.1
                )
smap.set_plot_params(cmap='Blues')
smap.set_shapefile()
for i, s in enumerate(sources):
    data = ods[s].copy().where(np.isfinite(ods[s]), 0)
    smap.set_data(data * 0)
    ax = grid[i]
    smap.set_topography(data)
    smap.visualize(ax=ax, addcbar=False, title=s)

# take care of uneven grids
if ax != grid[-1] and not grid[-1].title.get_text():
    grid[-1].remove()
    grid[-1].cax.remove()
if ax != grid[-2] and not grid[-2].title.get_text():
    grid[-2].remove()
    grid[-2].cax.remove()

plt.savefig(os.path.join(plot_dir, 'dem_topo_shade.png'), dpi=150, bbox_inches='tight')
../../_images/442f84e78ffb61a2bb5977c21d7cf47979a4d673f4bc63f6bcd23a23bcaece30.png

Slope#

fig = plt.figure(figsize=(x_size, y_size))
grid = AxesGrid(fig, 111,
                nrows_ncols=(n_rows, n_cols),
                axes_pad=0.7,
                cbar_mode='each',
                cbar_location='right',
                cbar_pad=0.1
                )

smap.set_topography()
smap.set_plot_params(vmin=0, vmax=0.7, cmap='Blues')

for i, s in enumerate(sources):
    data = ods[s + '_slope']
    smap.set_data(data)
    ax = grid[i]
    smap.visualize(ax=ax, addcbar=False, title=s + ' (slope)')
    cax = grid.cbar_axes[i]
    smap.colorbarbase(cax)

# take care of uneven grids
if ax != grid[-1] and not grid[-1].title.get_text():
    grid[-1].remove()
    grid[-1].cax.remove()
if ax != grid[-2] and not grid[-2].title.get_text():
    grid[-2].remove()
    grid[-2].cax.remove()

plt.savefig(os.path.join(plot_dir, 'dem_slope.png'), dpi=150, bbox_inches='tight'):
  Cell In[14], line 29
    plt.savefig(os.path.join(plot_dir, 'dem_slope.png'), dpi=150, bbox_inches='tight'):
                                                                                      ^
SyntaxError: invalid syntax

Some simple statistics about the DEMs#

df = pd.DataFrame()
for s in sources:
    df[s] = ods[s].data.flatten()[ods.mask.data.flatten() == 1]

dfs = pd.DataFrame()
for s in sources:
    dfs[s] = ods[s + '_slope'].data.flatten()[ods.mask.data.flatten() == 1]
dfs = df.describe()
dfs.loc['range'] = dfs.loc['max'] - dfs.loc['min']
dfs
ARCTICDEM COPDEM90 ALASKA DEM3 MAPZEN TANDEM COPDEM30 AW3D30 ASTER
count 767.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 0.0 1671.000000
mean 5529.413086 5277.694824 5316.735352 5328.549372 5319.226212 5153.806152 5277.514160 NaN 5295.587672
std 393.222046 370.293671 351.459656 351.724101 350.840381 388.357605 370.142792 NaN 351.229245
min 4871.422852 4584.852539 4610.906738 4644.000000 4622.000000 4497.636719 4585.657715 NaN 4565.000000
25% 5111.490234 4980.320312 5047.515381 5058.000000 5050.000000 4837.172607 4982.135010 NaN 5026.500000
50% 5613.539551 5231.637207 5255.337402 5259.000000 5258.000000 5044.266113 5231.779785 NaN 5250.000000
75% 5869.918945 5564.726562 5578.190430 5590.000000 5578.500000 5441.289551 5562.813232 NaN 5543.000000
max 6126.963379 6073.276855 6116.865234 6120.000000 6111.000000 5975.219727 6073.898438 NaN 6113.000000
range 1255.540527 1488.424316 1505.958496 1476.000000 1489.000000 1477.583008 1488.240723 NaN 1548.000000

Comparison matrix plot#

# Table of differences between DEMS
df_diff = pd.DataFrame()
done = []
for s1, s2 in itertools.product(sources, sources):
    if s1 == s2:
        continue
    if (s2, s1) in done:
        continue
    df_diff[s1 + '-' + s2] = df[s1] - df[s2]
    done.append((s1, s2))
# Decide on plot levels
max_diff = df_diff.quantile(0.99).max()
base_levels = np.array([-8, -5, -3, -1.5, -1, -0.5, -0.2, -0.1, 0, 0.1, 0.2, 0.5, 1, 1.5, 3, 5, 8])
if max_diff < 10:
    levels = base_levels
elif max_diff < 100:
    levels = base_levels * 10
elif max_diff < 1000:
    levels = base_levels * 100
else:
    levels = base_levels * 1000
levels = [l for l in levels if abs(l) < max_diff]
if max_diff > 10:
    levels = [int(l) for l in levels]
levels
[-800,
 -500,
 -300,
 -150,
 -100,
 -50,
 -20,
 -10,
 0,
 10,
 20,
 50,
 100,
 150,
 300,
 500,
 800]
smap.set_plot_params(levels=levels, cmap='PuOr', extend='both')
smap.set_shapefile(gdir.read_shapefile('outlines'))

fig = plt.figure(figsize=(14, 14))
grid = AxesGrid(fig, 111,
                nrows_ncols=(ns - 1, ns - 1),
                axes_pad=0.3,
                cbar_mode='single',
                cbar_location='right',
                cbar_pad=0.1
                )
done = []
for ax in grid:
    ax.set_axis_off()
for s1, s2 in itertools.product(sources, sources):
    if s1 == s2:
        continue
    if (s2, s1) in done:
        continue
    data = ods[s1] - ods[s2]
    ax = grid[sources.index(s1) * (ns - 1) + sources[1:].index(s2)]
    ax.set_axis_on()
    smap.set_data(data)
    smap.visualize(ax=ax, addcbar=False)
    done.append((s1, s2))
    ax.set_title(s1 + '-' + s2, fontsize=8)

cax = grid.cbar_axes[0]
smap.colorbarbase(cax);

plt.savefig(os.path.join(plot_dir, 'dem_diffs.png'), dpi=150, bbox_inches='tight')
../../_images/adf7bc68da82468a9bfd68fc8fbe697a04ad8c5cdd151ae8cbb3cdd66e11e93e.png

Comparison scatter plot#

import seaborn as sns
sns.set(style="ticks")

l1, l2 = (utils.nicenumber(df.min().min(), binsize=50, lower=True),
          utils.nicenumber(df.max().max(), binsize=50, lower=False))

def plot_unity():
    points = np.linspace(l1, l2, 100)
    plt.gca().plot(points, points, color='k', marker=None,
                   linestyle=':', linewidth=3.0)

g = sns.pairplot(df.dropna(how='all', axis=1).dropna(), plot_kws=dict(s=50, edgecolor="C0", linewidth=1));
g.map_offdiag(plot_unity)
for asx in g.axes:
    for ax in asx:
        ax.set_xlim((l1, l2))
        ax.set_ylim((l1, l2))

plt.savefig(os.path.join(plot_dir, 'dem_scatter.png'), dpi=150, bbox_inches='tight')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[20], line 13
      9     plt.gca().plot(points, points, color='k', marker=None,
     10                    linestyle=':', linewidth=3.0)
     11 
     12 g = sns.pairplot(df.dropna(how='all', axis=1).dropna(), plot_kws=dict(s=50, edgecolor="C0", linewidth=1));
---> 13 g.map_offdiag(plot_unity)
     14 for asx in g.axes:
     15     for ax in asx:
     16         ax.set_xlim((l1, l2))

File /usr/local/pyenv/versions/3.13.13/lib/python3.13/site-packages/seaborn/axisgrid.py:1425, in PairGrid.map_offdiag(self, func, **kwargs)
   1414 """Plot with a bivariate function on the off-diagonal subplots.
   1415 
   1416 Parameters
   (...)   1422 
   1423 """
   1424 if self.square_grid:
-> 1425     self.map_lower(func, **kwargs)
   1426     if not self._corner:
   1427         self.map_upper(func, **kwargs)

File /usr/local/pyenv/versions/3.13.13/lib/python3.13/site-packages/seaborn/axisgrid.py:1395, in PairGrid.map_lower(self, func, **kwargs)
   1384 """Plot with a bivariate function on the lower diagonal subplots.
   1385 
   1386 Parameters
   (...)   1392 
   1393 """
   1394 indices = zip(*np.tril_indices_from(self.axes, -1))
-> 1395 self._map_bivariate(func, indices, **kwargs)
   1396 return self

File /usr/local/pyenv/versions/3.13.13/lib/python3.13/site-packages/seaborn/axisgrid.py:1574, in PairGrid._map_bivariate(self, func, indices, **kwargs)
   1572     if ax is None:  # i.e. we are in corner mode
   1573         continue
-> 1574     self._plot_bivariate(x_var, y_var, ax, func, **kws)
   1575 self._add_axis_labels()
   1577 if "hue" in signature(func).parameters:

File /usr/local/pyenv/versions/3.13.13/lib/python3.13/site-packages/seaborn/axisgrid.py:1583, in PairGrid._plot_bivariate(self, x_var, y_var, ax, func, **kwargs)
   1581 """Draw a bivariate plot on the specified axes."""
   1582 if "hue" not in signature(func).parameters:
-> 1583     self._plot_bivariate_iter_hue(x_var, y_var, ax, func, **kwargs)
   1584     return
   1586 kwargs = kwargs.copy()

File /usr/local/pyenv/versions/3.13.13/lib/python3.13/site-packages/seaborn/axisgrid.py:1659, in PairGrid._plot_bivariate_iter_hue(self, x_var, y_var, ax, func, **kwargs)
   1657         func(x=x, y=y, **kws)
   1658     else:
-> 1659         func(x, y, **kws)
   1661 self._update_legend_data(ax)

TypeError: plot_unity() got an unexpected keyword argument 'color'
../../_images/e11aed5cfa4d037a1cc057179dd6efe43a86b5230c9d0fbd227a094928d2369c.png

Table statistics#

df.describe()
ARCTICDEM COPDEM90 ALASKA DEM3 MAPZEN TANDEM COPDEM30 AW3D30 ASTER
count 767.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 0.0 1671.000000
mean 5529.413086 5277.694824 5316.735352 5328.549372 5319.226212 5153.806152 5277.514160 NaN 5295.587672
std 393.222046 370.293671 351.459656 351.724101 350.840381 388.357605 370.142792 NaN 351.229245
min 4871.422852 4584.852539 4610.906738 4644.000000 4622.000000 4497.636719 4585.657715 NaN 4565.000000
25% 5111.490234 4980.320312 5047.515381 5058.000000 5050.000000 4837.172607 4982.135010 NaN 5026.500000
50% 5613.539551 5231.637207 5255.337402 5259.000000 5258.000000 5044.266113 5231.779785 NaN 5250.000000
75% 5869.918945 5564.726562 5578.190430 5590.000000 5578.500000 5441.289551 5562.813232 NaN 5543.000000
max 6126.963379 6073.276855 6116.865234 6120.000000 6111.000000 5975.219727 6073.898438 NaN 6113.000000
df.corr()
ARCTICDEM COPDEM90 ALASKA DEM3 MAPZEN TANDEM COPDEM30 AW3D30 ASTER
ARCTICDEM 1.000000 0.957821 0.957513 0.957046 0.957611 0.858754 0.957143 NaN 0.960869
COPDEM90 0.957821 1.000000 0.996075 0.995559 0.996033 0.928387 0.999891 NaN 0.995702
ALASKA 0.957513 0.996075 1.000000 0.998666 0.999929 0.930687 0.995697 NaN 0.997280
DEM3 0.957046 0.995559 0.998666 1.000000 0.998911 0.934220 0.995141 NaN 0.997277
MAPZEN 0.957611 0.996033 0.999929 0.998911 1.000000 0.930963 0.995645 NaN 0.997390
TANDEM 0.858754 0.928387 0.930687 0.934220 0.930963 1.000000 0.928303 NaN 0.923299
COPDEM30 0.957143 0.999891 0.995697 0.995141 0.995645 0.928303 1.000000 NaN 0.995283
AW3D30 NaN NaN NaN NaN NaN NaN NaN NaN NaN
ASTER 0.960869 0.995702 0.997280 0.997277 0.997390 0.923299 0.995283 NaN 1.000000
df_diff.describe()
ARCTICDEM-COPDEM90 ARCTICDEM-ALASKA ARCTICDEM-DEM3 ARCTICDEM-MAPZEN ARCTICDEM-TANDEM ARCTICDEM-COPDEM30 ARCTICDEM-AW3D30 ARCTICDEM-ASTER COPDEM90-ALASKA COPDEM90-DEM3 ... MAPZEN-TANDEM MAPZEN-COPDEM30 MAPZEN-AW3D30 MAPZEN-ASTER TANDEM-COPDEM30 TANDEM-AW3D30 TANDEM-ASTER COPDEM30-AW3D30 COPDEM30-ASTER AW3D30-ASTER
count 767.000000 767.000000 767.000000 767.000000 767.000000 767.000000 0.0 767.000000 1671.000000 1671.000000 ... 1671.000000 1671.000000 0.0 1671.000000 1671.000000 0.0 1671.000000 0.0 1671.000000 0.0
mean 110.184074 85.975136 76.316795 83.946521 251.653717 110.319916 NaN 102.185113 -39.040920 -50.854447 ... 165.419989 41.711812 NaN 23.638540 -123.708183 NaN -141.781450 NaN -18.073272 NaN
std 113.430321 113.401772 114.381777 113.273900 224.423676 114.278862 NaN 108.948818 37.098961 38.750864 ... 142.198023 38.776733 NaN 25.367381 144.721756 NaN 149.341491 NaN 39.801988 NaN
min -131.313965 -147.653809 -151.276367 -146.276367 -48.799805 -136.359863 NaN -121.276367 -130.026367 -163.009766 ... -16.507324 -38.260254 NaN -54.000000 -812.099121 NaN -829.932129 NaN -128.154297 NaN
25% 36.502441 5.880127 -1.839355 4.010010 111.032227 35.619141 NaN 28.157715 -65.058350 -74.026611 ... 71.849121 13.209717 NaN 6.000000 -178.104980 NaN -193.438232 NaN -42.359863 NaN
50% 81.454102 58.945312 38.955078 56.379883 170.058594 82.882324 NaN 71.489746 -33.313477 -39.517090 ... 123.718750 34.193848 NaN 22.000000 -95.392578 NaN -100.195801 NaN -13.353027 NaN
75% 158.221924 143.909668 139.544922 143.172607 358.547607 158.155029 NaN 171.116211 -10.340332 -23.536377 ... 214.740479 69.063965 NaN 40.000000 -18.154541 NaN -44.948730 NaN 8.449463 NaN
max 452.215332 401.538574 399.646973 396.646973 1085.924805 459.774902 NaN 428.646973 42.941895 17.004395 ... 851.888672 137.803711 NaN 130.000000 84.786133 NaN 85.662109 NaN 123.904297 NaN

8 rows × 36 columns

df_diff.abs().describe()
ARCTICDEM-COPDEM90 ARCTICDEM-ALASKA ARCTICDEM-DEM3 ARCTICDEM-MAPZEN ARCTICDEM-TANDEM ARCTICDEM-COPDEM30 ARCTICDEM-AW3D30 ARCTICDEM-ASTER COPDEM90-ALASKA COPDEM90-DEM3 ... MAPZEN-TANDEM MAPZEN-COPDEM30 MAPZEN-AW3D30 MAPZEN-ASTER TANDEM-COPDEM30 TANDEM-AW3D30 TANDEM-ASTER COPDEM30-AW3D30 COPDEM30-ASTER AW3D30-ASTER
count 767.000000 767.000000 767.000000 767.000000 767.000000 767.000000 0.0 767.000000 1671.000000 1671.000000 ... 1671.000000 1671.000000 0.0 1671.000000 1671.000000 0.0 1671.000000 0.0 1671.000000 0.0
mean 117.711647 105.086456 95.275925 103.447394 253.502304 118.163139 NaN 112.711556 42.339176 51.431971 ... 165.650691 44.739207 NaN 27.788151 135.184494 NaN 147.605664 NaN 32.801881 NaN
std 105.587395 95.936180 99.124644 95.769833 222.330780 106.137970 NaN 98.004051 33.283024 37.980538 ... 141.929044 35.238571 NaN 20.735849 134.057739 NaN 143.584187 NaN 28.886490 NaN
min 0.624512 0.155273 0.184570 0.289062 0.293457 0.650879 NaN 0.155762 0.000000 0.002441 ... 0.151367 0.117188 NaN 0.000000 0.002441 NaN 0.175293 NaN 0.005371 NaN
25% 41.026611 36.557861 21.433594 33.356445 111.032227 41.187012 NaN 41.791260 14.249756 23.536377 ... 71.849121 16.361816 NaN 12.000000 38.805664 NaN 49.170898 NaN 11.155762 NaN
50% 82.970215 65.577148 52.738770 64.257812 170.058594 84.257812 NaN 73.289062 33.496582 39.517090 ... 123.718750 34.610840 NaN 24.000000 95.392578 NaN 100.195801 NaN 23.466309 NaN
75% 158.221924 144.205811 141.093506 144.176514 358.547607 158.155029 NaN 171.116211 65.058350 74.026611 ... 214.740479 69.063965 NaN 40.000000 178.104980 NaN 193.438232 NaN 47.615723 NaN
max 452.215332 401.538574 399.646973 396.646973 1085.924805 459.774902 NaN 428.646973 130.026367 163.009766 ... 851.888672 137.803711 NaN 130.000000 812.099121 NaN 829.932129 NaN 128.154297 NaN

8 rows × 36 columns

What’s next?#