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-03-09 23:43:11: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2026-03-09 23:43:11: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2026-03-09 23:43:11: oggm.cfg: Multiprocessing: using all available processors (N=4)
2026-03-09 23:43:11: 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/data/gdirs/dems_v2/default (high resolution version: https://cluster.klima.uni-bremen.de/data/gdirs/dems_v1/highres)

# URL of the preprocessed GDirs
gdir_url = 'https://cluster.klima.uni-bremen.de/data/gdirs/dems_v2/default'
# 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-03-09 23:43:11: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2026-03-09 23:43:11: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

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: ['COPDEM90', 'TANDEM', 'COPDEM30', 'DEM3', 'ASTER', 'NASADEM', 'MAPZEN', 'SRTM', 'AW3D30']
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()

	COPDEM90	TANDEM	COPDEM30	DEM3	ASTER	NASADEM	MAPZEN	SRTM	AW3D30
count	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000
mean	3014.512695	3066.820068	3014.124512	3052.045681	3031.316346	3028.332505	3023.202610	3032.272219	3023.871970
std	259.186768	258.874115	259.403259	237.228101	250.883864	247.615467	254.144292	246.447678	253.597569
min	2416.773438	2469.287598	2413.545654	2490.000000	2428.000000	2430.000000	2413.000000	2450.000000	2417.000000
25%	2836.975037	2889.410828	2837.251770	2883.500000	2865.000000	2855.000000	2850.000000	2857.000000	2851.000000
50%	3044.824097	3097.242310	3044.196289	3068.500000	3058.000000	3054.000000	3051.000000	3060.000000	3052.500000
75%	3196.834839	3249.972717	3196.260803	3217.000000	3203.000000	3202.000000	3201.750000	3203.000000	3201.750000
max	3688.330078	3738.977051	3694.253174	3701.000000	3693.000000	3691.000000	3723.000000	3684.000000	3720.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

[-80, -50, -30, -15, -10, -5, -2, -1, 0, 1, 2, 5, 10, 15, 30, 50, 80]

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

Table statistics

df.describe()

	COPDEM90	TANDEM	COPDEM30	DEM3	ASTER	NASADEM	MAPZEN	SRTM	AW3D30
count	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000
mean	3014.512695	3066.820068	3014.124512	3052.045681	3031.316346	3028.332505	3023.202610	3032.272219	3023.871970
std	259.186768	258.874115	259.403259	237.228101	250.883864	247.615467	254.144292	246.447678	253.597569
min	2416.773438	2469.287598	2413.545654	2490.000000	2428.000000	2430.000000	2413.000000	2450.000000	2417.000000
25%	2836.975037	2889.410828	2837.251770	2883.500000	2865.000000	2855.000000	2850.000000	2857.000000	2851.000000
50%	3044.824097	3097.242310	3044.196289	3068.500000	3058.000000	3054.000000	3051.000000	3060.000000	3052.500000
75%	3196.834839	3249.972717	3196.260803	3217.000000	3203.000000	3202.000000	3201.750000	3203.000000	3201.750000
max	3688.330078	3738.977051	3694.253174	3701.000000	3693.000000	3691.000000	3723.000000	3684.000000	3720.000000

df.corr()

	COPDEM90	TANDEM	COPDEM30	DEM3	ASTER	NASADEM	MAPZEN	SRTM	AW3D30
COPDEM90	1.000000	0.999974	0.999958	0.998282	0.999365	0.999379	0.999700	0.998394	0.999660
TANDEM	0.999974	1.000000	0.999932	0.998220	0.999344	0.999342	0.999687	0.998315	0.999633
COPDEM30	0.999958	0.999932	1.000000	0.998258	0.999343	0.999358	0.999754	0.998275	0.999679
DEM3	0.998282	0.998220	0.998258	1.000000	0.998119	0.999294	0.998774	0.997864	0.998638
ASTER	0.999365	0.999344	0.999343	0.998119	1.000000	0.999131	0.999455	0.998518	0.999315
NASADEM	0.999379	0.999342	0.999358	0.999294	0.999131	1.000000	0.999643	0.998544	0.999553
MAPZEN	0.999700	0.999687	0.999754	0.998774	0.999455	0.999643	1.000000	0.998217	0.999819
SRTM	0.998394	0.998315	0.998275	0.997864	0.998518	0.998544	0.998217	1.000000	0.998216
AW3D30	0.999660	0.999633	0.999679	0.998638	0.999315	0.999553	0.999819	0.998216	1.000000

df_diff.describe()

	COPDEM90-TANDEM	COPDEM90-COPDEM30	COPDEM90-DEM3	COPDEM90-ASTER	COPDEM90-NASADEM	COPDEM90-MAPZEN	COPDEM90-SRTM	COPDEM90-AW3D30	TANDEM-COPDEM30	TANDEM-DEM3	...	ASTER-NASADEM	ASTER-MAPZEN	ASTER-SRTM	ASTER-AW3D30	NASADEM-MAPZEN	NASADEM-SRTM	NASADEM-AW3D30	MAPZEN-SRTM	MAPZEN-AW3D30	SRTM-AW3D30
count	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	...	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000
mean	-52.307072	0.388321	-37.532874	-16.803539	-13.819699	-8.689804	-17.759413	-9.359164	52.695396	14.774199	...	2.983841	8.113735	-0.955873	7.444375	5.129894	-3.939714	4.460534	-9.069608	-0.669360	8.400249
std	1.883556	2.387619	26.332881	12.309234	14.616837	8.057372	19.169901	8.710851	3.063715	26.213201	...	10.892745	8.953622	14.247553	9.720664	9.360127	13.380009	9.588316	16.809936	4.859039	16.555147
min	-79.299805	-15.591797	-124.903320	-52.883789	-69.138428	-44.972168	-73.730957	-41.972168	23.009766	-72.397461	...	-40.000000	-30.000000	-69.000000	-33.000000	-51.000000	-69.000000	-48.000000	-84.000000	-35.000000	-75.000000
25%	-52.462891	-0.533203	-49.661804	-24.004944	-18.632935	-11.912292	-29.375549	-14.219543	51.849060	2.716553	...	-3.000000	2.000000	-10.000000	1.000000	-1.000000	-10.000000	-2.000000	-20.750000	-3.000000	-1.000000
50%	-52.428711	0.146362	-28.553101	-14.929077	-9.529785	-7.221680	-17.123169	-9.151611	52.562256	23.847168	...	4.000000	8.000000	0.000000	7.000000	3.000000	-4.000000	3.000000	-10.000000	-1.000000	9.000000
75%	-52.394043	1.260254	-19.681152	-8.315002	-4.234375	-4.038086	-6.643860	-4.334778	53.654480	32.586426	...	9.000000	14.000000	7.000000	13.000000	10.000000	1.000000	8.000000	-1.000000	1.000000	20.000000
max	-23.296631	15.551025	10.634766	21.495605	27.849121	32.682129	53.235596	29.402100	82.028564	63.065674	...	55.000000	46.000000	59.000000	44.000000	40.000000	59.000000	38.000000	69.000000	28.000000	79.000000

8 rows × 36 columns

df_diff.abs().describe()

	COPDEM90-TANDEM	COPDEM90-COPDEM30	COPDEM90-DEM3	COPDEM90-ASTER	COPDEM90-NASADEM	COPDEM90-MAPZEN	COPDEM90-SRTM	COPDEM90-AW3D30	TANDEM-COPDEM30	TANDEM-DEM3	...	ASTER-NASADEM	ASTER-MAPZEN	ASTER-SRTM	ASTER-AW3D30	NASADEM-MAPZEN	NASADEM-SRTM	NASADEM-AW3D30	MAPZEN-SRTM	MAPZEN-AW3D30	SRTM-AW3D30
count	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	...	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000
mean	52.307072	1.539755	37.564394	17.410436	14.487043	9.364936	21.517830	10.511262	52.695396	26.697562	...	8.758856	9.781231	10.924798	9.678682	7.451212	10.136109	7.379117	15.292107	3.354257	14.691734
std	1.883556	1.865470	26.287883	11.434442	13.955490	7.261413	14.826293	7.278676	3.063715	13.874268	...	7.128573	7.093540	9.193391	7.498083	7.641868	9.580093	7.574302	11.442731	3.578254	11.346579
min	23.296631	0.000732	0.001709	0.000977	0.001953	0.000488	0.012451	0.000977	23.009766	0.003906	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	52.394043	0.342957	19.681152	8.656433	4.906738	4.394409	10.496765	5.184937	51.849060	16.986389	...	3.000000	4.000000	4.000000	4.000000	2.000000	3.000000	2.000000	6.000000	1.000000	5.000000
50%	52.428711	0.886719	28.553101	15.060913	9.772339	7.369751	18.754150	9.402832	52.562256	28.024414	...	7.000000	9.000000	9.000000	8.000000	4.000000	7.000000	5.000000	13.000000	2.000000	13.000000
75%	52.462891	2.049927	49.661804	24.004944	18.664856	12.140320	29.730896	14.364990	53.654480	35.479187	...	12.000000	14.000000	16.000000	14.000000	10.000000	15.000000	10.000000	23.000000	4.000000	22.000000
max	79.299805	15.591797	124.903320	52.883789	69.138428	44.972168	73.730957	41.972168	82.028564	72.397461	...	55.000000	46.000000	69.000000	44.000000	51.000000	69.000000	48.000000	84.000000	35.000000	79.000000

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents