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 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-03-10 19:39:01: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-10 19:39:01: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-10 19:39:01: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-03-10 19:39:01: 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]

2023-03-10 19:39:02: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2023-03-10 19:39:02: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2023-03-10 19:39:02: 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', 'ASTER', 'TANDEM', 'DEM3', 'MAPZEN', 'SRTM', 'NASADEM']
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	ASTER	TANDEM	DEM3	MAPZEN	SRTM	NASADEM
count	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000
mean	3023.871970	3031.316346	3066.820068	3061.494096	3023.202610	3032.272219	3028.332505
std	253.597569	250.883864	258.874115	241.860064	254.144292	246.447678	247.615467
min	2417.000000	2428.000000	2469.287598	2501.000000	2413.000000	2450.000000	2430.000000
25%	2851.000000	2865.000000	2889.410828	2886.000000	2850.000000	2857.000000	2855.000000
50%	3052.500000	3058.000000	3097.242310	3076.000000	3051.000000	3060.000000	3054.000000
75%	3201.750000	3203.000000	3249.972717	3233.750000	3201.750000	3203.000000	3202.000000
max	3720.000000	3693.000000	3738.977051	3706.000000	3723.000000	3684.000000	3691.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(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()

	AW3D30	ASTER	TANDEM	DEM3	MAPZEN	SRTM	NASADEM
count	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000	3218.000000
mean	3023.871970	3031.316346	3066.820068	3061.494096	3023.202610	3032.272219	3028.332505
std	253.597569	250.883864	258.874115	241.860064	254.144292	246.447678	247.615467
min	2417.000000	2428.000000	2469.287598	2501.000000	2413.000000	2450.000000	2430.000000
25%	2851.000000	2865.000000	2889.410828	2886.000000	2850.000000	2857.000000	2855.000000
50%	3052.500000	3058.000000	3097.242310	3076.000000	3051.000000	3060.000000	3054.000000
75%	3201.750000	3203.000000	3249.972717	3233.750000	3201.750000	3203.000000	3202.000000
max	3720.000000	3693.000000	3738.977051	3706.000000	3723.000000	3684.000000	3691.000000

df.corr()

	AW3D30	ASTER	TANDEM	DEM3	MAPZEN	SRTM	NASADEM
AW3D30	1.000000	0.999315	0.999633	0.996423	0.999819	0.998216	0.999553
ASTER	0.999315	1.000000	0.999344	0.995739	0.999455	0.998518	0.999131
TANDEM	0.999633	0.999344	1.000000	0.996175	0.999687	0.998315	0.999342
DEM3	0.996423	0.995739	0.996175	1.000000	0.996358	0.997038	0.997409
MAPZEN	0.999819	0.999455	0.999687	0.996358	1.000000	0.998217	0.999643
SRTM	0.998216	0.998518	0.998315	0.997038	0.998217	1.000000	0.998544
NASADEM	0.999553	0.999131	0.999342	0.997409	0.999643	0.998544	1.000000

df_diff.describe()

	AW3D30-ASTER	AW3D30-TANDEM	AW3D30-DEM3	AW3D30-MAPZEN	AW3D30-SRTM	AW3D30-NASADEM	ASTER-TANDEM	ASTER-DEM3	ASTER-MAPZEN	ASTER-SRTM	...	TANDEM-DEM3	TANDEM-MAPZEN	TANDEM-SRTM	TANDEM-NASADEM	DEM3-MAPZEN	DEM3-SRTM	DEM3-NASADEM	MAPZEN-SRTM	MAPZEN-NASADEM	SRTM-NASADEM
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	-7.444375	-42.947910	-37.622126	0.669360	-8.400249	-4.460534	-35.503534	-30.177750	8.113735	-0.955873	...	5.325784	43.617270	34.547661	38.487375	38.291485	29.221877	33.161591	-9.069608	-5.129894	3.939714
std	9.720664	8.722370	24.010492	4.859039	16.555147	9.588316	12.211383	24.464741	8.953622	14.247553	...	27.721822	7.972539	19.220421	14.527865	24.467068	19.342194	18.533220	16.809936	9.360127	13.380009
min	-44.000000	-82.822998	-121.000000	-28.000000	-79.000000	-38.000000	-79.663330	-112.000000	-30.000000	-69.000000	...	-79.397461	1.023926	-21.226807	-16.636475	-100.000000	-47.000000	-91.000000	-84.000000	-40.000000	-59.000000
25%	-13.000000	-47.949890	-49.000000	-1.000000	-20.000000	-8.000000	-43.959106	-42.000000	2.000000	-10.000000	...	-7.438965	40.532715	22.963257	33.745605	25.000000	22.000000	25.000000	-20.750000	-10.000000	-1.000000
50%	-7.000000	-43.253540	-36.000000	1.000000	-9.000000	-3.000000	-37.429810	-29.000000	8.000000	0.000000	...	6.978149	45.182861	35.271484	42.692627	38.000000	32.000000	34.000000	-10.000000	-3.000000	4.000000
75%	-1.000000	-38.227905	-25.000000	3.000000	1.000000	2.000000	-28.404724	-20.000000	14.000000	7.000000	...	21.998291	48.300476	45.801941	48.055481	50.000000	40.000000	43.000000	-1.000000	1.000000	10.000000
max	33.000000	-2.486084	109.000000	35.000000	75.000000	48.000000	0.389160	94.000000	46.000000	59.000000	...	131.650635	85.858154	106.757324	80.386475	123.000000	84.000000	92.000000	69.000000	51.000000	69.000000

8 rows × 21 columns

df_diff.abs().describe()

	AW3D30-ASTER	AW3D30-TANDEM	AW3D30-DEM3	AW3D30-MAPZEN	AW3D30-SRTM	AW3D30-NASADEM	ASTER-TANDEM	ASTER-DEM3	ASTER-MAPZEN	ASTER-SRTM	...	TANDEM-DEM3	TANDEM-MAPZEN	TANDEM-SRTM	TANDEM-NASADEM	DEM3-MAPZEN	DEM3-SRTM	DEM3-NASADEM	MAPZEN-SRTM	MAPZEN-NASADEM	SRTM-NASADEM
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.678682	42.947910	39.443754	3.354257	14.691734	7.379117	35.503776	33.658173	9.781231	10.924798	...	21.395372	43.617270	35.206708	38.762453	40.133623	32.021753	34.688626	15.292107	7.451212	10.136109
std	7.498083	8.722370	20.882402	3.578254	11.346579	7.574302	12.210679	19.397900	7.093540	9.193391	...	18.411083	7.972539	17.984406	13.776885	21.310154	14.233504	15.486716	11.442731	7.641868	9.580093
min	0.000000	2.486084	0.000000	0.000000	0.000000	0.000000	0.389160	0.000000	0.000000	0.000000	...	0.001709	1.023926	0.166016	0.005859	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	4.000000	38.227905	26.000000	1.000000	5.000000	2.000000	28.404724	21.000000	4.000000	4.000000	...	7.188965	40.532715	22.963257	33.745605	26.000000	24.000000	25.000000	6.000000	2.000000	3.000000
50%	8.000000	43.253540	36.000000	2.000000	13.000000	5.000000	37.429810	30.000000	9.000000	9.000000	...	16.839111	45.182861	35.271484	42.692627	38.000000	32.000000	34.000000	13.000000	4.000000	7.000000
75%	14.000000	47.949890	49.000000	4.000000	22.000000	10.000000	43.959106	43.000000	14.000000	16.000000	...	30.157959	48.300476	45.801941	48.055481	50.000000	41.000000	43.000000	23.000000	10.000000	15.000000
max	44.000000	82.822998	121.000000	35.000000	79.000000	48.000000	79.663330	112.000000	46.000000	69.000000	...	131.650635	85.858154	106.757324	80.386475	123.000000	84.000000	92.000000	84.000000	51.000000	69.000000