Compare different DEMs for individual glaciers

For most glaciers in the world there are several digital elevation models (DEM) which cover the respective glacier. In OGGM we have currently implemented 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. For more information, visit the rgitools documentation about DEMs.

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 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
2022-01-16 14:53:37: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2022-01-16 14:53:37: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2022-01-16 14:53:37: oggm.cfg: Multiprocessing: using all available processors (N=2)
2022-01-16 14:53:37: 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.4/rgitopo/

# URL of the preprocessed GDirs
gdir_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/rgitopo/'
# We use OGGM to download the data
gdir = init_glacier_directories([rgi_id], from_prepro_level=1, prepro_border=10, 
                                 prepro_rgi_version='62', prepro_base_url=gdir_url)[0]
2022-01-16 14:53:38: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2022-01-16 14:53:38: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2022-01-16 14:53:38: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/rgitopo/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.4/rgitopo/RGI62/b_010/L1/RGI60-11/RGI60-11.00.tar...
2022-01-16 14:53:42: oggm.utils: /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/rgitopo/RGI62/b_010/L1/RGI60-11/RGI60-11.00.tar verified successfully.

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', 'DEM3', 'MAPZEN', 'COPDEM', 'NASADEM', 'SRTM', 'ASTER', 'TANDEM']
Plotting directory: /__w/tutorials/tutorials/notebooks/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 xr.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 xr.open_rasterio(gdir.get_filepath('glacier_mask')) as ds:
    ods['mask'] = ds.sel(band=1).load()
/tmp/ipykernel_779/682396704.py:5: DeprecationWarning: open_rasterio is Deprecated in favor of rioxarray. For information about transitioning, see: https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html
  with xr.open_rasterio(demfile) as ds:
/tmp/ipykernel_779/682396704.py:5: DeprecationWarning: open_rasterio is Deprecated in favor of rioxarray. For information about transitioning, see: https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html
  with xr.open_rasterio(demfile) as ds:
/tmp/ipykernel_779/682396704.py:12: DeprecationWarning: open_rasterio is Deprecated in favor of rioxarray. For information about transitioning, see: https://corteva.github.io/rioxarray/stable/getting_started/getting_started.html
  with xr.open_rasterio(gdir.get_filepath('glacier_mask')) as ds:
# 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]:
    grid[-1].remove()
    grid[-1].cax.remove()

plt.savefig(os.path.join(plot_dir, 'dem_topo_color.png'), dpi=150, bbox_inches='tight')
../_images/dem_comparison_19_0.png

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_mode='none',
                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]:
    grid[-1].remove()
    grid[-1].cax.remove()

plt.savefig(os.path.join(plot_dir, 'dem_topo_shade.png'), dpi=150, bbox_inches='tight')
../_images/dem_comparison_21_0.png

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]:
    grid[-1].remove()
    grid[-1].cax.remove()

plt.savefig(os.path.join(plot_dir, 'dem_slope.png'), dpi=150, bbox_inches='tight')
../_images/dem_comparison_23_0.png

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 DEM3 MAPZEN COPDEM NASADEM SRTM ASTER TANDEM
count 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000
mean 3023.871970 3061.494096 3023.202610 3014.512451 3028.332505 3032.272219 3031.316346 3066.819092
std 253.597569 241.860064 254.144292 259.186768 247.615467 246.447678 250.883864 258.874115
min 2417.000000 2501.000000 2413.000000 2416.773438 2430.000000 2450.000000 2428.000000 2469.287598
25% 2851.000000 2886.000000 2850.000000 2836.975037 2855.000000 2857.000000 2865.000000 2889.410828
50% 3052.500000 3076.000000 3051.000000 3044.824097 3054.000000 3060.000000 3058.000000 3097.242310
75% 3201.750000 3233.750000 3201.750000 3196.834839 3202.000000 3203.000000 3203.000000 3249.972717
max 3720.000000 3706.000000 3723.000000 3688.330078 3691.000000 3684.000000 3693.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')
../_images/dem_comparison_30_0.png

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')
../_images/dem_comparison_32_0.png

Table statistics

df.describe()
AW3D30 DEM3 MAPZEN COPDEM NASADEM SRTM ASTER TANDEM
count 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000
mean 3023.871970 3061.494096 3023.202610 3014.512451 3028.332505 3032.272219 3031.316346 3066.819092
std 253.597569 241.860064 254.144292 259.186768 247.615467 246.447678 250.883864 258.874115
min 2417.000000 2501.000000 2413.000000 2416.773438 2430.000000 2450.000000 2428.000000 2469.287598
25% 2851.000000 2886.000000 2850.000000 2836.975037 2855.000000 2857.000000 2865.000000 2889.410828
50% 3052.500000 3076.000000 3051.000000 3044.824097 3054.000000 3060.000000 3058.000000 3097.242310
75% 3201.750000 3233.750000 3201.750000 3196.834839 3202.000000 3203.000000 3203.000000 3249.972717
max 3720.000000 3706.000000 3723.000000 3688.330078 3691.000000 3684.000000 3693.000000 3738.977051
df.corr()
AW3D30 DEM3 MAPZEN COPDEM NASADEM SRTM ASTER TANDEM
AW3D30 1.000000 0.996423 0.999819 0.999660 0.999553 0.998216 0.999315 0.999633
DEM3 0.996423 1.000000 0.996358 0.996189 0.997409 0.997038 0.995739 0.996175
MAPZEN 0.999819 0.996358 1.000000 0.999700 0.999643 0.998217 0.999455 0.999687
COPDEM 0.999660 0.996189 0.999700 1.000000 0.999379 0.998394 0.999365 0.999974
NASADEM 0.999553 0.997409 0.999643 0.999379 1.000000 0.998544 0.999131 0.999342
SRTM 0.998216 0.997038 0.998217 0.998394 0.998544 1.000000 0.998518 0.998315
ASTER 0.999315 0.995739 0.999455 0.999365 0.999131 0.998518 1.000000 0.999344
TANDEM 0.999633 0.996175 0.999687 0.999974 0.999342 0.998315 0.999344 1.000000
df_diff.describe()
AW3D30-DEM3 AW3D30-MAPZEN AW3D30-COPDEM AW3D30-NASADEM AW3D30-SRTM AW3D30-ASTER AW3D30-TANDEM DEM3-MAPZEN DEM3-COPDEM DEM3-NASADEM ... COPDEM-NASADEM COPDEM-SRTM COPDEM-ASTER COPDEM-TANDEM NASADEM-SRTM NASADEM-ASTER NASADEM-TANDEM SRTM-ASTER SRTM-TANDEM ASTER-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 -37.622126 0.669360 9.359164 -4.460534 -8.400249 -7.444375 -42.947910 38.291485 46.981290 33.161591 ... -13.819699 -17.759413 -16.803539 -52.307034 -3.939714 -2.983841 -38.487375 0.955873 -34.547661 -35.503534
std 24.010492 4.859039 8.710851 9.588316 16.555147 9.720664 8.722370 24.467068 27.893263 18.533220 ... 14.616837 19.169901 12.309234 1.883556 13.380009 10.892745 14.527865 14.247553 19.220421 12.211383
min -121.000000 -28.000000 -29.402100 -38.000000 -79.000000 -44.000000 -82.822998 -100.000000 -86.843262 -91.000000 ... -69.138428 -73.730957 -52.883789 -79.299805 -69.000000 -55.000000 -80.386475 -59.000000 -106.757324 -79.663330
25% -49.000000 -1.000000 4.334778 -8.000000 -20.000000 -13.000000 -47.949890 25.000000 30.414978 25.000000 ... -18.632935 -29.375549 -24.004944 -52.462891 -10.000000 -9.000000 -48.055481 -7.000000 -45.801941 -43.959106
50% -36.000000 1.000000 9.151611 -3.000000 -9.000000 -7.000000 -43.253540 38.000000 45.352905 34.000000 ... -9.529785 -17.123169 -14.929077 -52.428711 -4.000000 -4.000000 -42.692627 0.000000 -35.271484 -37.429810
75% -25.000000 3.000000 14.219543 2.000000 1.000000 -1.000000 -38.227905 50.000000 59.555908 43.000000 ... -4.234375 -6.643860 -8.315002 -52.394043 1.000000 3.000000 -33.745605 10.000000 -22.963257 -28.404724
max 109.000000 35.000000 41.972168 48.000000 75.000000 33.000000 -2.486084 123.000000 131.903320 92.000000 ... 27.849121 53.235596 21.495605 -23.296631 59.000000 40.000000 16.636475 69.000000 21.226807 0.389160

8 rows × 28 columns

df_diff.abs().describe()
AW3D30-DEM3 AW3D30-MAPZEN AW3D30-COPDEM AW3D30-NASADEM AW3D30-SRTM AW3D30-ASTER AW3D30-TANDEM DEM3-MAPZEN DEM3-COPDEM DEM3-NASADEM ... COPDEM-NASADEM COPDEM-SRTM COPDEM-ASTER COPDEM-TANDEM NASADEM-SRTM NASADEM-ASTER NASADEM-TANDEM SRTM-ASTER SRTM-TANDEM ASTER-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 39.443754 3.354257 10.511262 7.379117 14.691734 9.678682 42.947910 40.133623 48.232542 34.688626 ... 14.487043 21.517830 17.410436 52.307034 10.136109 8.758856 38.762453 10.924798 35.206708 35.503776
std 20.882402 3.578254 7.278676 7.574302 11.346579 7.498083 8.722370 21.310154 25.668278 15.486716 ... 13.955490 14.826293 11.434442 1.883556 9.580093 7.128573 13.776885 9.193391 17.984406 12.210679
min 0.000000 0.000000 0.000977 0.000000 0.000000 0.000000 2.486084 0.000000 0.087646 0.000000 ... 0.001953 0.012451 0.000977 23.296631 0.000000 0.000000 0.005859 0.000000 0.166016 0.389160
25% 26.000000 1.000000 5.184937 2.000000 5.000000 4.000000 38.227905 26.000000 30.705994 25.000000 ... 4.906738 10.496765 8.656433 52.394043 3.000000 3.000000 33.745605 4.000000 22.963257 28.404724
50% 36.000000 2.000000 9.402832 5.000000 13.000000 8.000000 43.253540 38.000000 45.544189 34.000000 ... 9.772339 18.754150 15.060913 52.428711 7.000000 7.000000 42.692627 9.000000 35.271484 37.429810
75% 49.000000 4.000000 14.364990 10.000000 22.000000 14.000000 47.949890 50.000000 59.782166 43.000000 ... 18.664856 29.730896 24.004944 52.462891 15.000000 12.000000 48.055481 16.000000 45.801941 43.959106
max 121.000000 35.000000 41.972168 48.000000 79.000000 44.000000 82.822998 123.000000 131.903320 92.000000 ... 69.138428 73.730957 52.883789 79.299805 69.000000 55.000000 80.386475 69.000000 106.757324 79.663330

8 rows × 28 columns

What’s next?