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-07-20 13:38:20: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2026-07-20 13:38:20: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2026-07-20 13:38:20: oggm.cfg: Multiprocessing: using all available processors (N=4)
2026-07-20 13:38:20: 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-07-20 13:38:21: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2026-07-20 13:38:21: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2026-07-20 13:38:21: 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.008.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.008.tar...
2026-07-20 13:38:21: 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: ['SRTM', 'NASADEM', 'MAPZEN', 'TANDEM', 'ASTER', 'DEM3', 'COPDEM30', 'AW3D30', 'COPDEM90']
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()

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

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

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

df.corr()

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

df_diff.describe()

	SRTM-NASADEM	SRTM-MAPZEN	SRTM-TANDEM	SRTM-ASTER	SRTM-DEM3	SRTM-COPDEM30	SRTM-AW3D30	SRTM-COPDEM90	NASADEM-MAPZEN	NASADEM-TANDEM	...	ASTER-DEM3	ASTER-COPDEM30	ASTER-AW3D30	ASTER-COPDEM90	DEM3-COPDEM30	DEM3-AW3D30	DEM3-COPDEM90	COPDEM30-AW3D30	COPDEM30-COPDEM90	AW3D30-COPDEM90
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	3.909599	9.074868	-34.515398	1.004349	-19.857409	18.168124	8.423734	17.773239	5.165269	-38.424998	...	-20.861758	17.163774	7.419385	16.768890	38.025533	28.281143	37.630648	-9.744389	-0.394885	9.349505
std	13.306268	16.760635	19.189930	14.211992	18.295258	19.636149	16.492592	19.143038	9.311976	14.470125	...	20.212678	12.606325	9.677552	12.362052	26.544060	20.759226	26.316103	8.719621	2.296632	8.700743
min	-63.000000	-72.000000	-107.228760	-59.000000	-74.000000	-63.182861	-73.000000	-53.959717	-41.000000	-75.664795	...	-106.000000	-32.396973	-29.000000	-32.612549	-9.135010	-24.000000	-4.963135	-40.435303	-10.730957	-29.469238
25%	-1.000000	0.000000	-45.617188	-7.000000	-32.000000	7.017456	-0.500000	6.779907	-1.000000	-47.947998	...	-27.000000	8.674438	1.000000	8.283325	20.157104	15.000000	19.852539	-14.411499	-1.276245	4.221680
50%	4.000000	10.000000	-35.234131	0.000000	-21.000000	17.448486	9.000000	17.155762	3.000000	-42.638916	...	-18.000000	15.394531	7.000000	14.818359	29.390381	22.000000	28.607910	-9.586670	-0.166992	9.158936
75%	10.000000	21.000000	-22.934814	10.000000	-9.000000	30.036133	20.000000	29.323608	10.000000	-33.795532	...	-10.000000	24.747314	13.000000	24.107666	49.742188	36.000000	49.807861	-4.713745	0.518311	14.035889
max	71.000000	84.000000	21.623047	68.000000	58.000000	80.028076	82.000000	74.127686	39.000000	15.975098	...	46.000000	52.590332	46.000000	52.431641	131.028076	114.000000	125.127686	26.607422	16.592773	41.233398

8 rows × 36 columns

df_diff.abs().describe()

	SRTM-NASADEM	SRTM-MAPZEN	SRTM-TANDEM	SRTM-ASTER	SRTM-DEM3	SRTM-COPDEM30	SRTM-AW3D30	SRTM-COPDEM90	NASADEM-MAPZEN	NASADEM-TANDEM	...	ASTER-DEM3	ASTER-COPDEM30	ASTER-AW3D30	ASTER-COPDEM90	DEM3-COPDEM30	DEM3-AW3D30	DEM3-COPDEM90	COPDEM30-AW3D30	COPDEM30-COPDEM90	AW3D30-COPDEM90
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	10.039453	15.263747	35.187296	10.919851	23.011805	22.003013	14.641193	21.490858	7.433054	38.691837	...	22.702081	17.911928	9.676608	17.410026	38.091287	28.555763	37.664915	10.836687	1.499186	10.491179
std	9.566846	11.412421	17.927794	9.149292	14.122605	15.214420	11.338245	14.847274	7.624547	13.740463	...	18.120593	11.518214	7.419808	11.441161	26.449586	20.379699	26.267020	7.317387	1.783883	7.283151
min	0.000000	0.000000	0.001465	0.000000	0.000000	0.000244	0.000000	0.002197	0.000000	0.058594	...	0.000000	0.014893	0.000000	0.004639	0.022705	0.000000	0.172363	0.002197	0.001465	0.001709
25%	3.000000	6.000000	22.934814	4.000000	12.000000	10.488037	5.000000	10.327148	2.000000	33.795532	...	11.000000	9.050293	4.000000	8.570312	20.157104	15.000000	19.852539	5.573120	0.335815	5.152466
50%	7.000000	13.000000	35.234131	9.000000	22.000000	19.133545	12.000000	18.616943	5.000000	42.638916	...	19.000000	15.547363	8.000000	14.985840	29.390381	22.000000	28.607910	9.769287	0.867920	9.456299
75%	15.000000	23.000000	45.617188	16.000000	32.000000	30.468506	22.000000	29.641357	10.000000	47.947998	...	27.000000	24.809326	14.000000	24.158325	49.742188	36.000000	49.807861	14.518677	1.980347	14.206177
max	71.000000	84.000000	107.228760	68.000000	74.000000	80.028076	82.000000	74.127686	41.000000	75.664795	...	106.000000	52.590332	46.000000	52.431641	131.028076	114.000000	125.127686	40.435303	16.592773	41.233398

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents