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

2023-08-28 00:05:22: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-08-28 00:05:22: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-08-28 00:05:22: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-08-28 00:05:23: 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]

2023-08-28 00:05:24: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2023-08-28 00:05:24: 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: ['AW3D30', 'COPDEM90', 'DEM3', 'ASTER', 'MAPZEN', 'SRTM', 'COPDEM30', 'NASADEM', 'TANDEM']
Plotting directory: /__w/tutorials/tutorials/notebooks/others/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()

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

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

/usr/local/pyenv/versions/3.10.12/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

Table statistics

df.describe()

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

df.corr()

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

df_diff.describe()

	AW3D30-COPDEM90	AW3D30-DEM3	AW3D30-ASTER	AW3D30-MAPZEN	AW3D30-SRTM	AW3D30-COPDEM30	AW3D30-NASADEM	AW3D30-TANDEM	COPDEM90-DEM3	COPDEM90-ASTER	...	MAPZEN-SRTM	MAPZEN-COPDEM30	MAPZEN-NASADEM	MAPZEN-TANDEM	SRTM-COPDEM30	SRTM-NASADEM	SRTM-TANDEM	COPDEM30-NASADEM	COPDEM30-TANDEM	NASADEM-TANDEM
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	9.359164	-28.173710	-7.444375	0.669360	-8.400249	9.747485	-4.460534	-42.947910	-37.532874	-16.803539	...	-9.069608	9.078125	-5.129894	-43.617270	18.147733	3.939714	-34.547661	-14.208019	-52.695396	-38.487375
std	8.710851	20.781199	9.720664	4.859039	16.555147	8.717534	9.588316	8.722370	26.332881	12.309234	...	16.809936	7.750302	9.360127	7.972539	19.706512	13.380009	19.220421	14.878318	3.063715	14.527865
min	-29.402100	-114.000000	-44.000000	-28.000000	-79.000000	-23.341553	-38.000000	-82.822998	-124.903320	-52.883789	...	-84.000000	-31.172363	-40.000000	-85.858154	-55.908691	-59.000000	-106.757324	-71.462646	-82.028564	-80.386475
25%	4.334778	-36.000000	-13.000000	-1.000000	-20.000000	4.652039	-8.000000	-47.949890	-49.661804	-24.004944	...	-20.750000	4.448425	-10.000000	-48.300476	6.655701	-1.000000	-45.801941	-19.370178	-53.654480	-48.055481
50%	9.151611	-22.000000	-7.000000	1.000000	-9.000000	9.645752	-3.000000	-43.253540	-28.553101	-14.929077	...	-10.000000	7.459595	-3.000000	-45.182861	17.336670	4.000000	-35.271484	-10.123413	-52.562256	-42.692627
75%	14.219543	-15.000000	-1.000000	3.000000	1.000000	14.405151	2.000000	-38.227905	-19.681152	-8.315002	...	-1.000000	11.907959	1.000000	-40.532715	30.145203	10.000000	-22.963257	-4.427490	-51.849060	-33.745605
max	41.972168	27.000000	33.000000	35.000000	75.000000	43.652588	48.000000	-2.486084	10.634766	21.495605	...	69.000000	46.652588	51.000000	-1.023926	79.334717	69.000000	21.226807	27.088867	-23.009766	16.636475

8 rows × 36 columns

df_diff.abs().describe()

	AW3D30-COPDEM90	AW3D30-DEM3	AW3D30-ASTER	AW3D30-MAPZEN	AW3D30-SRTM	AW3D30-COPDEM30	AW3D30-NASADEM	AW3D30-TANDEM	COPDEM90-DEM3	COPDEM90-ASTER	...	MAPZEN-SRTM	MAPZEN-COPDEM30	MAPZEN-NASADEM	MAPZEN-TANDEM	SRTM-COPDEM30	SRTM-NASADEM	SRTM-TANDEM	COPDEM30-NASADEM	COPDEM30-TANDEM	NASADEM-TANDEM
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	10.511262	28.465817	9.678682	3.354257	14.691734	10.837928	7.379117	42.947910	37.564394	17.410436	...	15.292107	9.462017	7.451212	43.617270	22.048134	10.136109	35.206708	14.978851	52.695396	38.762453
std	7.278676	20.379114	7.498083	3.578254	11.346579	7.317186	7.574302	8.722370	26.287883	11.434442	...	11.442731	7.276508	7.641868	7.972539	15.215713	9.580093	17.984406	14.101744	3.063715	13.776885
min	0.000977	0.000000	0.000000	0.000000	0.000000	0.003418	0.000000	2.486084	0.001709	0.000977	...	0.000000	0.007568	0.000000	1.023926	0.003906	0.000000	0.166016	0.002197	23.009766	0.005859
25%	5.184937	15.000000	4.000000	1.000000	5.000000	5.617004	2.000000	38.227905	19.681152	8.656433	...	6.000000	4.625854	2.000000	40.532715	10.525879	3.000000	22.963257	5.027039	51.849060	33.745605
50%	9.402832	22.000000	8.000000	2.000000	13.000000	9.821289	5.000000	43.253540	28.553101	15.060913	...	13.000000	7.555664	4.000000	45.182861	19.279541	7.000000	35.271484	10.341919	52.562256	42.692627
75%	14.364990	36.000000	14.000000	4.000000	22.000000	14.540771	10.000000	47.949890	49.661804	24.004944	...	23.000000	12.061157	10.000000	48.300476	30.630066	15.000000	45.801941	19.433960	53.654480	48.055481
max	41.972168	114.000000	44.000000	35.000000	79.000000	43.652588	48.000000	82.822998	124.903320	52.883789	...	84.000000	46.652588	51.000000	85.858154	79.334717	69.000000	106.757324	71.462646	82.028564	80.386475

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents