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
/tmp/ipykernel_25941/1578510794.py:1: DeprecationWarning: 
Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),
(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)
but was not found to be installed on your system.
If this would cause problems for you,
please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466
        
  import pandas as pd
# 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-02-02 17:21:08: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-02-02 17:21:08: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-02-02 17:21:08: oggm.cfg: Multiprocessing: using all available processors (N=4)
2024-02-02 17:21:08: 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-02-02 17:21:08: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2024-02-02 17:21:08: 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', 'ASTER', 'COPDEM30', 'NASADEM', 'SRTM', 'AW3D30', 'DEM3', 'TANDEM', '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')
../../_images/f1dab1bd4d6a8acc47b69f23886054afbd4c881172d7bb82a5fb9e37a6d66d7d.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_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')
../../_images/335075c3ea3f720eacebe83368743d526b3d6f59fbec6d995e3c969e0db25155.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] 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')
../../_images/12c0e993a8ed4e2833bd0bcc08402b0f086b25f1dbd10807ef5c1c12e08f2d37.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()
COPDEM90 ASTER COPDEM30 NASADEM SRTM AW3D30 DEM3 TANDEM MAPZEN
count 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000
mean 3014.512695 3031.316346 3014.124512 3028.332505 3032.272219 3023.871970 3052.045681 3066.820068 3023.202610
std 259.186768 250.883864 259.403259 247.615467 246.447678 253.597569 237.228101 258.874115 254.144292
min 2416.773438 2428.000000 2413.545654 2430.000000 2450.000000 2417.000000 2490.000000 2469.287598 2413.000000
25% 2836.975037 2865.000000 2837.251770 2855.000000 2857.000000 2851.000000 2883.500000 2889.410828 2850.000000
50% 3044.824097 3058.000000 3044.196289 3054.000000 3060.000000 3052.500000 3068.500000 3097.242310 3051.000000
75% 3196.834839 3203.000000 3196.260803 3202.000000 3203.000000 3201.750000 3217.000000 3249.972717 3201.750000
max 3688.330078 3693.000000 3694.253174 3691.000000 3684.000000 3720.000000 3701.000000 3738.977051 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
[-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')
../../_images/cc7c73d336d33ef1174c30fad827667d9a376f60bb8d9bb478dc98f58617018b.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/18ea4dcb9d4e8c77d244a8a169e83557eee804f145ba938d1f0278d411d2e465.png

Table statistics#

df.describe()
COPDEM90 ASTER COPDEM30 NASADEM SRTM AW3D30 DEM3 TANDEM MAPZEN
count 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000 3218.000000
mean 3014.512695 3031.316346 3014.124512 3028.332505 3032.272219 3023.871970 3052.045681 3066.820068 3023.202610
std 259.186768 250.883864 259.403259 247.615467 246.447678 253.597569 237.228101 258.874115 254.144292
min 2416.773438 2428.000000 2413.545654 2430.000000 2450.000000 2417.000000 2490.000000 2469.287598 2413.000000
25% 2836.975037 2865.000000 2837.251770 2855.000000 2857.000000 2851.000000 2883.500000 2889.410828 2850.000000
50% 3044.824097 3058.000000 3044.196289 3054.000000 3060.000000 3052.500000 3068.500000 3097.242310 3051.000000
75% 3196.834839 3203.000000 3196.260803 3202.000000 3203.000000 3201.750000 3217.000000 3249.972717 3201.750000
max 3688.330078 3693.000000 3694.253174 3691.000000 3684.000000 3720.000000 3701.000000 3738.977051 3723.000000
df.corr()
COPDEM90 ASTER COPDEM30 NASADEM SRTM AW3D30 DEM3 TANDEM MAPZEN
COPDEM90 1.000000 0.999365 0.999958 0.999379 0.998394 0.999660 0.998282 0.999974 0.999700
ASTER 0.999365 1.000000 0.999343 0.999131 0.998518 0.999315 0.998119 0.999344 0.999455
COPDEM30 0.999958 0.999343 1.000000 0.999358 0.998275 0.999679 0.998258 0.999932 0.999754
NASADEM 0.999379 0.999131 0.999358 1.000000 0.998544 0.999553 0.999294 0.999342 0.999643
SRTM 0.998394 0.998518 0.998275 0.998544 1.000000 0.998216 0.997864 0.998315 0.998217
AW3D30 0.999660 0.999315 0.999679 0.999553 0.998216 1.000000 0.998638 0.999633 0.999819
DEM3 0.998282 0.998119 0.998258 0.999294 0.997864 0.998638 1.000000 0.998220 0.998774
TANDEM 0.999974 0.999344 0.999932 0.999342 0.998315 0.999633 0.998220 1.000000 0.999687
MAPZEN 0.999700 0.999455 0.999754 0.999643 0.998217 0.999819 0.998774 0.999687 1.000000
df_diff.describe()
COPDEM90-ASTER COPDEM90-COPDEM30 COPDEM90-NASADEM COPDEM90-SRTM COPDEM90-AW3D30 COPDEM90-DEM3 COPDEM90-TANDEM COPDEM90-MAPZEN ASTER-COPDEM30 ASTER-NASADEM ... SRTM-AW3D30 SRTM-DEM3 SRTM-TANDEM SRTM-MAPZEN AW3D30-DEM3 AW3D30-TANDEM AW3D30-MAPZEN DEM3-TANDEM DEM3-MAPZEN TANDEM-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 -16.803539 0.388321 -13.819699 -17.759413 -9.359164 -37.532874 -52.307072 -8.689804 17.191860 2.983841 ... 8.400249 -19.773462 -34.547661 9.069608 -28.173710 -42.947910 0.669360 -14.774199 28.843070 43.617270
std 12.309234 2.387619 14.616837 19.169901 8.710851 26.332881 1.883556 8.057372 12.572292 10.892745 ... 16.555147 18.297883 19.220421 16.809936 20.781199 8.722370 4.859039 26.213201 20.831143 7.972539
min -52.883789 -15.591797 -69.138428 -73.730957 -41.972168 -124.903320 -79.299805 -44.972168 -23.341553 -40.000000 ... -75.000000 -74.000000 -106.757324 -69.000000 -114.000000 -82.822998 -28.000000 -63.065674 -22.000000 1.023926
25% -24.004944 -0.533203 -18.632935 -29.375549 -14.219543 -49.661804 -52.462891 -11.912292 8.794128 -3.000000 ... -1.000000 -32.000000 -45.801941 1.000000 -36.000000 -47.949890 -1.000000 -32.586426 15.000000 40.532715
50% -14.929077 0.146362 -9.529785 -17.123169 -9.151611 -28.553101 -52.428711 -7.221680 15.391968 4.000000 ... 9.000000 -21.000000 -35.271484 10.000000 -22.000000 -43.253540 1.000000 -23.847168 23.000000 45.182861
75% -8.315002 1.260254 -4.234375 -6.643860 -4.334778 -19.681152 -52.394043 -4.038086 24.869507 9.000000 ... 20.000000 -9.000000 -22.963257 20.750000 -15.000000 -38.227905 3.000000 -2.716553 37.000000 48.300476
max 21.495605 15.551025 27.849121 53.235596 29.402100 10.634766 -23.296631 32.682129 52.308594 55.000000 ... 79.000000 54.000000 21.226807 84.000000 27.000000 -2.486084 35.000000 72.397461 116.000000 85.858154

8 rows × 36 columns

df_diff.abs().describe()
COPDEM90-ASTER COPDEM90-COPDEM30 COPDEM90-NASADEM COPDEM90-SRTM COPDEM90-AW3D30 COPDEM90-DEM3 COPDEM90-TANDEM COPDEM90-MAPZEN ASTER-COPDEM30 ASTER-NASADEM ... SRTM-AW3D30 SRTM-DEM3 SRTM-TANDEM SRTM-MAPZEN AW3D30-DEM3 AW3D30-TANDEM AW3D30-MAPZEN DEM3-TANDEM DEM3-MAPZEN TANDEM-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 17.410436 1.539755 14.487043 21.517830 10.511262 37.564394 52.307072 9.364936 17.918605 8.758856 ... 14.691734 22.972343 35.206708 15.292107 28.465817 42.947910 3.354257 26.697562 29.093536 43.617270
std 11.434442 1.865470 13.955490 14.826293 7.278676 26.287883 1.883556 7.261413 11.512526 7.128573 ... 11.346579 14.072358 17.984406 11.442731 20.379114 8.722370 3.578254 13.874268 20.479767 7.972539
min 0.000977 0.000732 0.001953 0.012451 0.000977 0.001709 23.296631 0.000488 0.000732 0.000000 ... 0.000000 0.000000 0.166016 0.000000 0.000000 2.486084 0.000000 0.003906 0.000000 1.023926
25% 8.656433 0.342957 4.906738 10.496765 5.184937 19.681152 52.394043 4.394409 9.204651 3.000000 ... 5.000000 12.000000 22.963257 6.000000 15.000000 38.227905 1.000000 16.986389 15.000000 40.532715
50% 15.060913 0.886719 9.772339 18.754150 9.402832 28.553101 52.428711 7.369751 15.604980 7.000000 ... 13.000000 22.000000 35.271484 13.000000 22.000000 43.253540 2.000000 28.024414 23.000000 45.182861
75% 24.004944 2.049927 18.664856 29.730896 14.364990 49.661804 52.462891 12.140320 24.869507 12.000000 ... 22.000000 32.000000 45.801941 23.000000 36.000000 47.949890 4.000000 35.479187 37.000000 48.300476
max 52.883789 15.591797 69.138428 73.730957 41.972168 124.903320 79.299805 44.972168 52.308594 55.000000 ... 79.000000 74.000000 106.757324 84.000000 114.000000 82.822998 35.000000 72.397461 116.000000 85.858154

8 rows × 36 columns

What’s next?#