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

2024-08-25 21:35:05: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:35:05: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:35:05: oggm.cfg: Multiprocessing: using all available processors (N=4)
2024-08-25 21:35:07: 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]

2024-08-25 21:35:07: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2024-08-25 21:35:07: 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: ['COPDEM30', 'TANDEM', 'DEM3', 'AW3D30', 'ASTER', 'COPDEM90', 'NASADEM', 'SRTM', 'MAPZEN']
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()

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

[-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(xdata, ydata, **kwargs):
    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')

Table statistics

df.describe()

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

df.corr()

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

df_diff.describe()

	COPDEM30-TANDEM	COPDEM30-DEM3	COPDEM30-AW3D30	COPDEM30-ASTER	COPDEM30-COPDEM90	COPDEM30-NASADEM	COPDEM30-SRTM	COPDEM30-MAPZEN	TANDEM-DEM3	TANDEM-AW3D30	...	ASTER-COPDEM90	ASTER-NASADEM	ASTER-SRTM	ASTER-MAPZEN	COPDEM90-NASADEM	COPDEM90-SRTM	COPDEM90-MAPZEN	NASADEM-SRTM	NASADEM-MAPZEN	SRTM-MAPZEN
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.695396	-37.921195	-9.747485	-17.191860	-0.388321	-14.208019	-18.147733	-9.078125	14.774199	42.947910	...	16.803539	2.983841	-0.955873	8.113735	-13.819699	-17.759413	-8.689804	-3.939714	5.129894	9.069608
std	3.063715	26.573851	8.717534	12.572292	2.387619	14.878318	19.706512	7.750302	26.213201	8.722370	...	12.309234	10.892745	14.247553	8.953622	14.616837	19.169901	8.057372	13.380009	9.360127	16.809936
min	-82.028564	-130.334717	-43.652588	-52.308594	-15.551025	-71.462646	-79.334717	-46.652588	-72.397461	2.486084	...	-21.495605	-40.000000	-69.000000	-30.000000	-69.138428	-73.730957	-44.972168	-69.000000	-51.000000	-69.000000
25%	-53.654480	-49.488220	-14.405151	-24.869507	-1.260254	-19.370178	-30.145203	-11.907959	2.716553	38.227905	...	8.315002	-3.000000	-10.000000	2.000000	-18.632935	-29.375549	-11.912292	-10.000000	-1.000000	1.000000
50%	-52.562256	-29.365234	-9.645752	-15.391968	-0.146362	-10.123413	-17.336670	-7.459595	23.847168	43.253540	...	14.929077	4.000000	0.000000	8.000000	-9.529785	-17.123169	-7.221680	-4.000000	3.000000	10.000000
75%	-51.849060	-19.947693	-4.652039	-8.794128	0.533203	-4.427490	-6.655701	-4.448425	32.586426	47.949890	...	24.004944	9.000000	7.000000	14.000000	-4.234375	-6.643860	-4.038086	1.000000	10.000000	20.750000
max	-23.009766	11.646729	23.341553	23.341553	15.591797	27.088867	55.908691	31.172363	63.065674	82.822998	...	52.883789	55.000000	59.000000	46.000000	27.849121	53.235596	32.682129	59.000000	40.000000	84.000000

8 rows × 36 columns

df_diff.abs().describe()

	COPDEM30-TANDEM	COPDEM30-DEM3	COPDEM30-AW3D30	COPDEM30-ASTER	COPDEM30-COPDEM90	COPDEM30-NASADEM	COPDEM30-SRTM	COPDEM30-MAPZEN	TANDEM-DEM3	TANDEM-AW3D30	...	ASTER-COPDEM90	ASTER-NASADEM	ASTER-SRTM	ASTER-MAPZEN	COPDEM90-NASADEM	COPDEM90-SRTM	COPDEM90-MAPZEN	NASADEM-SRTM	NASADEM-MAPZEN	SRTM-MAPZEN
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.695396	37.997709	10.837928	17.918605	1.539755	14.978851	22.048134	9.462017	26.697562	42.947910	...	17.410436	8.758856	10.924798	9.781231	14.487043	21.517830	9.364936	10.136109	7.451212	15.292107
std	3.063715	26.464294	7.317186	11.512526	1.865470	14.101744	15.215713	7.276508	13.874268	8.722370	...	11.434442	7.128573	9.193391	7.093540	13.955490	14.826293	7.261413	9.580093	7.641868	11.442731
min	23.009766	0.135986	0.003418	0.000732	0.000732	0.002197	0.003906	0.007568	0.003906	2.486084	...	0.000977	0.000000	0.000000	0.000000	0.001953	0.012451	0.000488	0.000000	0.000000	0.000000
25%	51.849060	19.947693	5.617004	9.204651	0.342957	5.027039	10.525879	4.625854	16.986389	38.227905	...	8.656433	3.000000	4.000000	4.000000	4.906738	10.496765	4.394409	3.000000	2.000000	6.000000
50%	52.562256	29.365234	9.821289	15.604980	0.886719	10.341919	19.279541	7.555664	28.024414	43.253540	...	15.060913	7.000000	9.000000	9.000000	9.772339	18.754150	7.369751	7.000000	4.000000	13.000000
75%	53.654480	49.488220	14.540771	24.869507	2.049927	19.433960	30.630066	12.061157	35.479187	47.949890	...	24.004944	12.000000	16.000000	14.000000	18.664856	29.730896	12.140320	15.000000	10.000000	23.000000
max	82.028564	130.334717	43.652588	52.308594	15.591797	71.462646	79.334717	46.652588	72.397461	82.822998	...	52.883789	55.000000	69.000000	46.000000	69.138428	73.730957	44.972168	69.000000	51.000000	84.000000

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents