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-03-09 23:43:33: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2026-03-09 23:43:33: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2026-03-09 23:43:33: oggm.cfg: Multiprocessing: using all available processors (N=4)
2026-03-09 23:43:33: 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/2025.4'
# 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-03-09 23:43:34: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2026-03-09 23:43:34: 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: ['COPDEM90', 'ARCTICDEM', 'TANDEM', 'COPDEM30', 'DEM3', 'ASTER', 'MAPZEN', 'ALASKA', 'AW3D30']
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/28cf901f9c36a69a5774e31cb57520f35faa2f778f8cfa09fb17b631302a9c98.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/5229419f12aa8dea143e7e7c784eca50952c3085f488cd44f81fcd7123449369.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
COPDEM90 ARCTICDEM TANDEM COPDEM30 DEM3 ASTER MAPZEN ALASKA AW3D30
count 1671.000000 767.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 0.0
mean 5277.694824 5529.413086 5153.806152 5277.514160 5328.549372 5295.587672 5319.226212 5316.735352 NaN
std 370.293671 393.222046 388.357605 370.142792 351.724101 351.229245 350.840381 351.459656 NaN
min 4584.852539 4871.422852 4497.636719 4585.657715 4644.000000 4565.000000 4622.000000 4610.906738 NaN
25% 4980.320312 5111.490234 4837.172607 4982.135010 5058.000000 5026.500000 5050.000000 5047.515381 NaN
50% 5231.637207 5613.539551 5044.266113 5231.779785 5259.000000 5250.000000 5258.000000 5255.337402 NaN
75% 5564.726562 5869.918945 5441.289551 5562.813232 5590.000000 5543.000000 5578.500000 5578.190430 NaN
max 6073.276855 6126.963379 5975.219727 6073.898438 6120.000000 6113.000000 6111.000000 6116.865234 NaN
range 1488.424316 1255.540527 1477.583008 1488.240723 1476.000000 1548.000000 1489.000000 1505.958496 NaN

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/51fbad2a315edb38ee8ed28893cb02fa3acd004821aa3d085b7e5b4bb363292f.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)
     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:

File /usr/local/pyenv/versions/3.13.12/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.12/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.12/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.12/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.12/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/a8d2bb06bdc11193e8e6759a41467312ce26a057b1a3356371ea08f8f2378a56.png

Table statistics#

df.describe()
COPDEM90 ARCTICDEM TANDEM COPDEM30 DEM3 ASTER MAPZEN ALASKA AW3D30
count 1671.000000 767.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 0.0
mean 5277.694824 5529.413086 5153.806152 5277.514160 5328.549372 5295.587672 5319.226212 5316.735352 NaN
std 370.293671 393.222046 388.357605 370.142792 351.724101 351.229245 350.840381 351.459656 NaN
min 4584.852539 4871.422852 4497.636719 4585.657715 4644.000000 4565.000000 4622.000000 4610.906738 NaN
25% 4980.320312 5111.490234 4837.172607 4982.135010 5058.000000 5026.500000 5050.000000 5047.515381 NaN
50% 5231.637207 5613.539551 5044.266113 5231.779785 5259.000000 5250.000000 5258.000000 5255.337402 NaN
75% 5564.726562 5869.918945 5441.289551 5562.813232 5590.000000 5543.000000 5578.500000 5578.190430 NaN
max 6073.276855 6126.963379 5975.219727 6073.898438 6120.000000 6113.000000 6111.000000 6116.865234 NaN
df.corr()
COPDEM90 ARCTICDEM TANDEM COPDEM30 DEM3 ASTER MAPZEN ALASKA AW3D30
COPDEM90 1.000000 0.957821 0.928387 0.999891 0.995559 0.995702 0.996033 0.996075 NaN
ARCTICDEM 0.957821 1.000000 0.858754 0.957143 0.957046 0.960869 0.957611 0.957513 NaN
TANDEM 0.928387 0.858754 1.000000 0.928303 0.934220 0.923299 0.930963 0.930687 NaN
COPDEM30 0.999891 0.957143 0.928303 1.000000 0.995141 0.995283 0.995645 0.995697 NaN
DEM3 0.995559 0.957046 0.934220 0.995141 1.000000 0.997277 0.998911 0.998666 NaN
ASTER 0.995702 0.960869 0.923299 0.995283 0.997277 1.000000 0.997390 0.997280 NaN
MAPZEN 0.996033 0.957611 0.930963 0.995645 0.998911 0.997390 1.000000 0.999929 NaN
ALASKA 0.996075 0.957513 0.930687 0.995697 0.998666 0.997280 0.999929 1.000000 NaN
AW3D30 NaN NaN NaN NaN NaN NaN NaN NaN NaN
df_diff.describe()
COPDEM90-ARCTICDEM COPDEM90-TANDEM COPDEM90-COPDEM30 COPDEM90-DEM3 COPDEM90-ASTER COPDEM90-MAPZEN COPDEM90-ALASKA COPDEM90-AW3D30 ARCTICDEM-TANDEM ARCTICDEM-COPDEM30 ... DEM3-ASTER DEM3-MAPZEN DEM3-ALASKA DEM3-AW3D30 ASTER-MAPZEN ASTER-ALASKA ASTER-AW3D30 MAPZEN-ALASKA MAPZEN-AW3D30 ALASKA-AW3D30
count 767.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 0.0 767.000000 767.000000 ... 1671.000000 1671.000000 1671.000000 0.0 1671.000000 1671.000000 0.0 1671.000000 0.0 0.0
mean -110.184074 123.888710 0.180525 -50.854447 -17.892747 -41.531287 -39.040920 NaN 251.653717 110.319916 ... 32.961700 9.323160 11.813529 NaN -23.638540 -21.148171 NaN 2.490369 NaN NaN
std 113.430321 144.647827 5.457573 38.750864 38.488241 37.537576 37.098961 NaN 224.423676 114.278862 ... 25.941984 16.420294 18.159255 NaN 25.367381 25.915689 NaN 4.234587 NaN NaN
min -452.215332 -78.665039 -17.572266 -163.009766 -124.032227 -135.499023 -130.026367 NaN -48.799805 -136.359863 ... -52.000000 -47.000000 -44.825195 NaN -130.000000 -137.751953 NaN -14.705078 NaN NaN
25% -158.221924 17.190430 -2.751221 -74.026611 -40.641357 -68.542236 -65.058350 NaN 111.032227 35.619141 ... 19.000000 0.000000 1.023926 NaN -40.000000 -36.852539 NaN 0.054199 NaN NaN
50% -81.454102 95.284180 0.013672 -39.517090 -12.738770 -34.726074 -33.313477 NaN 170.058594 82.882324 ... 29.000000 10.000000 11.001953 NaN -22.000000 -19.398438 NaN 2.204590 NaN NaN
75% -36.502441 179.401855 3.210938 -23.536377 7.423828 -12.118164 -10.340332 NaN 358.547607 158.155029 ... 48.000000 20.000000 23.316650 NaN -6.000000 -1.874023 NaN 5.209473 NaN NaN
max 131.313965 812.890137 18.569824 17.004395 111.338867 33.757324 42.941895 NaN 1085.924805 459.774902 ... 128.000000 60.000000 73.005371 NaN 54.000000 55.279785 NaN 20.434082 NaN NaN

8 rows × 36 columns

df_diff.abs().describe()
COPDEM90-ARCTICDEM COPDEM90-TANDEM COPDEM90-COPDEM30 COPDEM90-DEM3 COPDEM90-ASTER COPDEM90-MAPZEN COPDEM90-ALASKA COPDEM90-AW3D30 ARCTICDEM-TANDEM ARCTICDEM-COPDEM30 ... DEM3-ASTER DEM3-MAPZEN DEM3-ALASKA DEM3-AW3D30 ASTER-MAPZEN ASTER-ALASKA ASTER-AW3D30 MAPZEN-ALASKA MAPZEN-AW3D30 ALASKA-AW3D30
count 767.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 1671.000000 0.0 767.000000 767.000000 ... 1671.000000 1671.000000 1671.000000 0.0 1671.000000 1671.000000 0.0 1671.000000 0.0 0.0
mean 117.711647 134.962799 4.060298 51.431971 31.869296 44.107472 42.339176 NaN 253.502304 118.163139 ... 35.226212 15.171155 17.200776 NaN 27.788151 26.224301 NaN 3.799071 NaN NaN
std 105.587395 134.367798 3.649907 37.980538 28.025463 34.471963 33.283024 NaN 222.330780 106.137970 ... 22.770194 11.238222 13.166538 NaN 20.735849 20.760783 NaN 3.113807 NaN NaN
min 0.624512 0.050781 0.000488 0.002441 0.000000 0.069336 0.000000 NaN 0.293457 0.650879 ... 0.000000 0.000000 0.002441 NaN 0.000000 0.001953 NaN 0.002930 NaN NaN
25% 41.026611 38.032227 1.359131 23.536377 10.204102 14.985596 14.249756 NaN 111.032227 41.187012 ... 20.000000 6.000000 6.530273 NaN 12.000000 11.052246 NaN 1.328125 NaN NaN
50% 82.970215 95.284180 2.911621 39.517090 22.945801 34.726074 33.496582 NaN 170.058594 84.257812 ... 29.000000 12.000000 13.911133 NaN 24.000000 21.200195 NaN 3.072266 NaN NaN
75% 158.221924 179.401855 5.748779 74.026611 46.310059 68.542236 65.058350 NaN 358.547607 158.155029 ... 48.000000 22.000000 25.122559 NaN 40.000000 37.241211 NaN 5.534912 NaN NaN
max 452.215332 812.890137 18.569824 163.009766 124.032227 135.499023 130.026367 NaN 1085.924805 459.774902 ... 128.000000 60.000000 73.005371 NaN 130.000000 137.751953 NaN 20.434082 NaN NaN

8 rows × 36 columns

What’s next?#