Compare different DEMs for individual glaciers: RGI-TOPO for RGI v6.0

For most glaciers in the world there are several digital elevation models (DEM) which cover the respective glacier. In OGGM we have currently implemented more than 10 different open access DEMs to choose from. Some are regional and only available in certain areas (e.g. Greenland or Antarctica) and some cover almost the entire globe.

This notebook allows to see which of the DEMs are available for a selected glacier and how they compare to each other. That way it is easy to spot systematic differences and also invalid points in the DEMs.

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 = 'RGI60-11.00897'

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

import pandas as pd
import numpy as np
from oggm import cfg, utils, workflow, tasks, graphics, GlacierDirectory
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-DEMS', reset=True)
cfg.PARAMS['use_intersects'] = False

2026-05-15 17:18:05: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2026-05-15 17:18:05: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2026-05-15 17:18:05: oggm.cfg: Multiprocessing: using all available processors (N=4)
2026-05-15 17:18:05: 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 https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2026.1/all_dems/ (version with the selected DEM: https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2026.1/selected_dem/)

# 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_base_url=gdir_url)[0]

2026-05-15 17:18:05: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2026-05-15 17:18:06: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2026-05-15 17:18:06: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2026.1/all_dems/RGI62/b_010/L1/RGI60-11/RGI60-11.00.tar to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2026.1/all_dems/RGI62/b_010/L1/RGI60-11/RGI60-11.00.tar...

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: RGI60-11.00897
Available DEM sources: ['TANDEM', 'MAPZEN', 'SRTM', 'AW3D30', 'NASADEM', 'ASTER', 'DEM3', 'COPDEM90', 'COPDEM30']
Plotting directory: /__w/tutorials/tutorials/notebooks/tutorials/outputs/RGI60-11.00897

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

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

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

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]

df.describe()

	TANDEM	MAPZEN	SRTM	AW3D30	NASADEM	ASTER	DEM3	COPDEM90	COPDEM30
count	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000
mean	3066.144775	3022.554520	3031.629388	3023.205654	3027.719789	3030.625039	3051.486797	3013.856201	3013.461426
std	259.755615	255.076764	247.360252	254.510913	248.589552	251.752913	238.175459	260.099152	260.305847
min	2469.645020	2412.000000	2450.000000	2416.000000	2431.000000	2428.000000	2488.000000	2417.131104	2413.272217
25%	2887.448853	2848.500000	2856.000000	2850.000000	2853.500000	2864.000000	2882.000000	2834.996582	2835.018066
50%	3095.948730	3051.000000	3059.000000	3052.000000	3052.000000	3056.000000	3067.000000	3043.559570	3042.966797
75%	3248.545532	3200.500000	3204.500000	3203.000000	3203.000000	3205.000000	3217.000000	3196.595459	3195.463745
max	3737.720215	3716.000000	3684.000000	3719.000000	3692.000000	3691.000000	3703.000000	3688.516113	3694.617676

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

[-100, -50, -20, -10, 0, 10, 20, 50, 100]

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

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,  marker=None,
                   linestyle=':', linewidth=3.0, color='k')

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[19], line 13
      9     plt.gca().plot(points, points,  marker=None,
     10                    linestyle=':', linewidth=3.0, color='k')
     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'

Table statistics

df.describe()

	TANDEM	MAPZEN	SRTM	AW3D30	NASADEM	ASTER	DEM3	COPDEM90	COPDEM30
count	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000
mean	3066.144775	3022.554520	3031.629388	3023.205654	3027.719789	3030.625039	3051.486797	3013.856201	3013.461426
std	259.755615	255.076764	247.360252	254.510913	248.589552	251.752913	238.175459	260.099152	260.305847
min	2469.645020	2412.000000	2450.000000	2416.000000	2431.000000	2428.000000	2488.000000	2417.131104	2413.272217
25%	2887.448853	2848.500000	2856.000000	2850.000000	2853.500000	2864.000000	2882.000000	2834.996582	2835.018066
50%	3095.948730	3051.000000	3059.000000	3052.000000	3052.000000	3056.000000	3067.000000	3043.559570	3042.966797
75%	3248.545532	3200.500000	3204.500000	3203.000000	3203.000000	3205.000000	3217.000000	3196.595459	3195.463745
max	3737.720215	3716.000000	3684.000000	3719.000000	3692.000000	3691.000000	3703.000000	3688.516113	3694.617676

df.corr()

	TANDEM	MAPZEN	SRTM	AW3D30	NASADEM	ASTER	DEM3	COPDEM90	COPDEM30
TANDEM	1.000000	0.999687	0.998330	0.999636	0.999344	0.999343	0.998222	0.999974	0.999935
MAPZEN	0.999687	1.000000	0.998246	0.999822	0.999648	0.999463	0.998783	0.999701	0.999751
SRTM	0.998330	0.998246	1.000000	0.998246	0.998573	0.998533	0.997875	0.998413	0.998307
AW3D30	0.999636	0.999822	0.998246	1.000000	0.999564	0.999329	0.998646	0.999664	0.999680
NASADEM	0.999344	0.999648	0.998573	0.999564	1.000000	0.999137	0.999301	0.999384	0.999368
ASTER	0.999343	0.999463	0.998533	0.999329	0.999137	1.000000	0.998130	0.999365	0.999346
DEM3	0.998222	0.998783	0.997875	0.998646	0.999301	0.998130	1.000000	0.998290	0.998267
COPDEM90	0.999974	0.999701	0.998413	0.999664	0.999384	0.999365	0.998290	1.000000	0.999961
COPDEM30	0.999935	0.999751	0.998307	0.999680	0.999368	0.999346	0.998267	0.999961	1.000000

df_diff.describe()

	TANDEM-MAPZEN	TANDEM-SRTM	TANDEM-AW3D30	TANDEM-NASADEM	TANDEM-ASTER	TANDEM-DEM3	TANDEM-COPDEM90	TANDEM-COPDEM30	MAPZEN-SRTM	MAPZEN-AW3D30	...	NASADEM-ASTER	NASADEM-DEM3	NASADEM-COPDEM90	NASADEM-COPDEM30	ASTER-DEM3	ASTER-COPDEM90	ASTER-COPDEM30	DEM3-COPDEM90	DEM3-COPDEM30	COPDEM90-COPDEM30
count	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	...	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000
mean	43.590266	34.515398	42.939132	38.424998	35.519748	14.657989	52.288639	52.683521	-9.074868	-0.651134	...	-2.905250	-23.767008	13.863640	14.258524	-20.861758	16.768890	17.163774	37.630648	38.025533	0.394885
std	7.957814	19.189930	8.697157	14.470125	12.246588	26.184944	1.918336	3.006184	16.760635	4.841786	...	10.862782	13.827482	14.565806	14.798913	20.212678	12.362052	12.606325	26.316103	26.544060	2.296632
min	-0.142090	-21.623047	-3.142090	-15.975098	0.063477	-72.623047	20.051270	17.101074	-84.000000	-31.000000	...	-51.000000	-85.000000	-18.129395	-18.503906	-106.000000	-32.612549	-32.396973	-4.963135	-9.135010	-16.592773
25%	40.418091	22.934814	38.299927	33.795532	28.249023	2.590820	52.394165	51.883301	-21.000000	-3.000000	...	-9.000000	-31.000000	4.275879	4.520020	-27.000000	8.283325	8.674438	19.852539	20.157104	-0.518311
50%	45.169434	35.234131	43.225586	42.638916	37.508057	23.722900	52.427734	52.589600	-10.000000	-1.000000	...	-3.000000	-21.000000	9.619873	10.160156	-18.000000	14.818359	15.394531	28.607910	29.390381	0.166992
75%	48.326782	45.617188	48.114868	47.947998	44.028198	32.422241	52.462891	53.701782	0.000000	1.000000	...	3.000000	-15.000000	18.536743	19.309448	-10.000000	24.107666	24.747314	49.807861	49.742188	1.276245
max	83.664795	107.228760	83.232178	75.664795	74.897705	58.887939	76.923584	79.928711	72.000000	34.000000	...	40.000000	12.000000	68.478271	70.410156	46.000000	52.431641	52.590332	125.127686	131.028076	10.730957

8 rows × 36 columns

df_diff.abs().describe()

	TANDEM-MAPZEN	TANDEM-SRTM	TANDEM-AW3D30	TANDEM-NASADEM	TANDEM-ASTER	TANDEM-DEM3	TANDEM-COPDEM90	TANDEM-COPDEM30	MAPZEN-SRTM	MAPZEN-AW3D30	...	NASADEM-ASTER	NASADEM-DEM3	NASADEM-COPDEM90	NASADEM-COPDEM30	ASTER-DEM3	ASTER-COPDEM90	ASTER-COPDEM30	DEM3-COPDEM90	DEM3-COPDEM30	COPDEM90-COPDEM30
count	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	...	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000	3219.000000
mean	43.590355	35.187296	42.941478	38.691837	35.519748	26.621369	52.288639	52.683521	15.263747	3.359428	...	8.727555	23.908046	14.493817	14.973715	22.702081	17.410026	17.911928	37.664915	38.091287	1.499186
std	7.957330	17.927794	8.685566	13.740463	12.246588	13.844030	1.918336	3.006184	11.412421	3.546510	...	7.088670	13.582091	13.938685	14.074612	18.120593	11.441161	11.518214	26.267020	26.449586	1.783883
min	0.142090	0.001465	0.632812	0.058594	0.063477	0.021484	20.051270	17.101074	0.000000	0.000000	...	0.000000	0.000000	0.040283	0.003418	0.000000	0.004639	0.014893	0.172363	0.022705	0.001465
25%	40.418091	22.934814	38.299927	33.795532	28.249023	16.813110	52.394165	51.883301	6.000000	1.000000	...	3.000000	15.000000	4.831055	5.109253	11.000000	8.570312	9.050293	19.852539	20.157104	0.335815
50%	45.169434	35.234131	43.225586	42.638916	37.508057	27.969238	52.427734	52.589600	13.000000	2.000000	...	7.000000	21.000000	9.785400	10.431396	19.000000	14.985840	15.547363	28.607910	29.390381	0.867920
75%	48.326782	45.617188	48.114868	47.947998	44.028198	35.264526	52.462891	53.701782	23.000000	4.000000	...	12.000000	31.000000	18.536743	19.309448	27.000000	24.158325	24.809326	49.807861	49.742188	1.980347
max	83.664795	107.228760	83.232178	75.664795	74.897705	72.623047	76.923584	79.928711	84.000000	34.000000	...	51.000000	85.000000	68.478271	70.410156	106.000000	52.431641	52.590332	125.127686	131.028076	16.592773

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents