RGI-TOPO for RGI 7.0

OGGM was used to generate the topography data used to compute the topographical attributes and the centerlines products for RGI v7.0.

Here we show how to access this data from OGGM.

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 = 'RGI2000-v7.0-G-01-06486'  # Denali

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

from oggm import cfg, utils, workflow, tasks, graphics, GlacierDirectory
import pandas as pd
import numpy as np
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-RGITOPO-RGI7', reset=True)
cfg.PARAMS['use_intersects'] = False

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

# URL of the preprocessed GDirs
gdir_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2023.1/'
# We use OGGM to download the data
gdir = init_glacier_directories([rgi_id], from_prepro_level=1, prepro_border=10,  prepro_rgi_version='70', prepro_base_url=gdir_url)[0]

2025-06-05 13:08:22: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2025-06-05 13:08:22: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2025-06-05 13:08:22: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2023.1/RGI70/b_010/L1/RGI2000-v7.0-G-01/RGI2000-v7.0-G-01-06.tar to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/rgitopo/2023.1/RGI70/b_010/L1/RGI2000-v7.0-G-01/RGI2000-v7.0-G-01-06.tar...
2025-06-05 13:08:31: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~fmaussion/misc/rgi7_data/00_rgi70_regions/00_rgi70_O1Regions/00_rgi70_O1Regions.dbf to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~fmaussion/misc/rgi7_data/00_rgi70_regions/00_rgi70_O1Regions/00_rgi70_O1Regions.dbf...
2025-06-05 13:08:31: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~fmaussion/misc/rgi7_data/00_rgi70_regions/00_rgi70_O2Regions/00_rgi70_O2Regions.dbf to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~fmaussion/misc/rgi7_data/00_rgi70_regions/00_rgi70_O2Regions/00_rgi70_O2Regions.dbf...

gdir

<oggm.GlacierDirectory>
  RGI id: RGI2000-v7.0-G-01-06486
  Region: 01: Alaska
  Subregion: 01-02: Alaska Range (Wrangell/Kilbuck)
  Glacier type: Glacier
  Terminus type: Not assigned
  Status: Glacier
  Area: 0.961122833004822 km2
  Lon, Lat: (-151.0094740399913, 63.061062)
  Grid (nx, ny): (70, 88)
  Grid (dx, dy): (24.0, -24.0)

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: RGI2000-v7.0-G-01-06486
Available DEM sources: ['MAPZEN', 'ASTER', 'ARCTICDEM', 'DEM3', 'TANDEM', 'ALASKA', 'COPDEM90', 'AW3D30', 'COPDEM30']
Plotting directory: /__w/tutorials/tutorials/notebooks/tutorials/outputs/RGI2000-v7.0-G-01-06486

# 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]

dfs = df.describe()
dfs.loc['range'] = dfs.loc['max'] - dfs.loc['min']
dfs

	MAPZEN	ASTER	ARCTICDEM	DEM3	TANDEM	ALASKA	COPDEM90	AW3D30	COPDEM30
count	1671.000000	1671.000000	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	0.0	1671.000000
mean	5319.226212	5295.587672	5529.413086	5328.549372	5153.806152	5316.735352	5277.694824	NaN	5277.514160
std	350.840381	351.229245	393.222046	351.724101	388.357605	351.459656	370.293671	NaN	370.142792
min	4622.000000	4565.000000	4871.422852	4644.000000	4497.636719	4610.906738	4584.852539	NaN	4585.657715
25%	5050.000000	5026.500000	5111.490234	5058.000000	4837.172607	5047.515381	4980.320312	NaN	4982.135010
50%	5258.000000	5250.000000	5613.539551	5259.000000	5044.266113	5255.337402	5231.637207	NaN	5231.779785
75%	5578.500000	5543.000000	5869.918945	5590.000000	5441.289551	5578.190430	5564.726562	NaN	5562.813232
max	6111.000000	6113.000000	6126.963379	6120.000000	5975.219727	6116.865234	6073.276855	NaN	6073.898438
range	1489.000000	1548.000000	1255.540527	1476.000000	1477.583008	1505.958496	1488.424316	NaN	1488.240723

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

[-800,
 -500,
 -300,
 -150,
 -100,
 -50,
 -20,
 -10,
 0,
 10,
 20,
 50,
 100,
 150,
 300,
 500,
 800]

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

	MAPZEN	ASTER	ARCTICDEM	DEM3	TANDEM	ALASKA	COPDEM90	AW3D30	COPDEM30
count	1671.000000	1671.000000	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	0.0	1671.000000
mean	5319.226212	5295.587672	5529.413086	5328.549372	5153.806152	5316.735352	5277.694824	NaN	5277.514160
std	350.840381	351.229245	393.222046	351.724101	388.357605	351.459656	370.293671	NaN	370.142792
min	4622.000000	4565.000000	4871.422852	4644.000000	4497.636719	4610.906738	4584.852539	NaN	4585.657715
25%	5050.000000	5026.500000	5111.490234	5058.000000	4837.172607	5047.515381	4980.320312	NaN	4982.135010
50%	5258.000000	5250.000000	5613.539551	5259.000000	5044.266113	5255.337402	5231.637207	NaN	5231.779785
75%	5578.500000	5543.000000	5869.918945	5590.000000	5441.289551	5578.190430	5564.726562	NaN	5562.813232
max	6111.000000	6113.000000	6126.963379	6120.000000	5975.219727	6116.865234	6073.276855	NaN	6073.898438

df.corr()

	MAPZEN	ASTER	ARCTICDEM	DEM3	TANDEM	ALASKA	COPDEM90	AW3D30	COPDEM30
MAPZEN	1.000000	0.997390	0.957611	0.998911	0.930963	0.999929	0.996033	NaN	0.995645
ASTER	0.997390	1.000000	0.960869	0.997277	0.923299	0.997280	0.995702	NaN	0.995283
ARCTICDEM	0.957611	0.960869	1.000000	0.957046	0.858754	0.957513	0.957821	NaN	0.957143
DEM3	0.998911	0.997277	0.957046	1.000000	0.934220	0.998666	0.995559	NaN	0.995141
TANDEM	0.930963	0.923299	0.858754	0.934220	1.000000	0.930687	0.928387	NaN	0.928303
ALASKA	0.999929	0.997280	0.957513	0.998666	0.930687	1.000000	0.996075	NaN	0.995697
COPDEM90	0.996033	0.995702	0.957821	0.995559	0.928387	0.996075	1.000000	NaN	0.999891
AW3D30	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
COPDEM30	0.995645	0.995283	0.957143	0.995141	0.928303	0.995697	0.999891	NaN	1.000000

df_diff.describe()

	MAPZEN-ASTER	MAPZEN-ARCTICDEM	MAPZEN-DEM3	MAPZEN-TANDEM	MAPZEN-ALASKA	MAPZEN-COPDEM90	MAPZEN-AW3D30	MAPZEN-COPDEM30	ASTER-ARCTICDEM	ASTER-DEM3	...	TANDEM-ALASKA	TANDEM-COPDEM90	TANDEM-AW3D30	TANDEM-COPDEM30	ALASKA-COPDEM90	ALASKA-AW3D30	ALASKA-COPDEM30	COPDEM90-AW3D30	COPDEM90-COPDEM30	AW3D30-COPDEM30
count	1671.000000	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	0.0	1671.000000	767.000000	1671.000000	...	1671.000000	1671.000000	0.0	1671.000000	1671.000000	0.0	1671.000000	0.0	1671.000000	0.0
mean	23.638540	-83.946521	-9.323160	165.419989	2.490369	41.531287	NaN	41.711812	-102.185113	-32.961700	...	-162.929626	-123.888710	NaN	-123.708183	39.040920	NaN	39.221443	NaN	0.180525	NaN
std	25.367381	113.273900	16.420294	142.198023	4.234587	37.537576	NaN	38.776733	108.948818	25.941984	...	142.417313	144.647827	NaN	144.721756	37.098961	NaN	38.322636	NaN	5.457573	NaN
min	-54.000000	-396.646973	-60.000000	-16.507324	-14.705078	-33.757324	NaN	-38.260254	-428.646973	-128.000000	...	-848.312500	-812.890137	NaN	-812.099121	-42.941895	NaN	-46.030762	NaN	-17.572266	NaN
25%	6.000000	-143.172607	-20.000000	71.849121	0.054199	12.118164	NaN	13.209717	-171.116211	-48.000000	...	-212.847412	-179.401855	NaN	-178.104980	10.340332	NaN	10.385498	NaN	-2.751221	NaN
50%	22.000000	-56.379883	-10.000000	123.718750	2.204590	34.726074	NaN	34.193848	-71.489746	-29.000000	...	-120.674805	-95.284180	NaN	-95.392578	33.313477	NaN	32.329102	NaN	0.013672	NaN
75%	40.000000	-4.010010	0.000000	214.740479	5.209473	68.542236	NaN	69.063965	-28.157715	-19.000000	...	-67.788818	-17.190430	NaN	-18.154541	65.058350	NaN	65.968750	NaN	3.210938	NaN
max	130.000000	146.276367	47.000000	851.888672	20.434082	135.499023	NaN	137.803711	121.276367	52.000000	...	15.732910	78.665039	NaN	84.786133	130.026367	NaN	135.432617	NaN	18.569824	NaN

8 rows × 36 columns

df_diff.abs().describe()

	MAPZEN-ASTER	MAPZEN-ARCTICDEM	MAPZEN-DEM3	MAPZEN-TANDEM	MAPZEN-ALASKA	MAPZEN-COPDEM90	MAPZEN-AW3D30	MAPZEN-COPDEM30	ASTER-ARCTICDEM	ASTER-DEM3	...	TANDEM-ALASKA	TANDEM-COPDEM90	TANDEM-AW3D30	TANDEM-COPDEM30	ALASKA-COPDEM90	ALASKA-AW3D30	ALASKA-COPDEM30	COPDEM90-AW3D30	COPDEM90-COPDEM30	AW3D30-COPDEM30
count	1671.000000	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	0.0	1671.000000	767.000000	1671.000000	...	1671.000000	1671.000000	0.0	1671.000000	1671.000000	0.0	1671.000000	0.0	1671.000000	0.0
mean	27.788151	103.447394	15.171155	165.650691	3.799071	44.107472	NaN	44.739207	112.711556	35.226212	...	163.309280	134.962799	NaN	135.184494	42.339176	NaN	42.954536	NaN	4.060298	NaN
std	20.735849	95.769833	11.238222	141.929044	3.113807	34.471963	NaN	35.238571	98.004051	22.770194	...	141.981552	134.367798	NaN	134.057739	33.283024	NaN	34.083324	NaN	3.649907	NaN
min	0.000000	0.289062	0.000000	0.151367	0.002930	0.069336	NaN	0.117188	0.155762	0.000000	...	0.070801	0.050781	NaN	0.002441	0.000000	NaN	0.007324	NaN	0.000488	NaN
25%	12.000000	33.356445	6.000000	71.849121	1.328125	14.985596	NaN	16.361816	41.791260	20.000000	...	67.788818	38.032227	NaN	38.805664	14.249756	NaN	14.802979	NaN	1.359131	NaN
50%	24.000000	64.257812	12.000000	123.718750	3.072266	34.726074	NaN	34.610840	73.289062	29.000000	...	120.674805	95.284180	NaN	95.392578	33.496582	NaN	32.889648	NaN	2.911621	NaN
75%	40.000000	144.176514	22.000000	214.740479	5.534912	68.542236	NaN	69.063965	171.116211	48.000000	...	212.847412	179.401855	NaN	178.104980	65.058350	NaN	65.968750	NaN	5.748779	NaN
max	130.000000	396.646973	60.000000	851.888672	20.434082	135.499023	NaN	137.803711	428.646973	128.000000	...	848.312500	812.890137	NaN	812.099121	130.026367	NaN	135.432617	NaN	18.569824	NaN

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents