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

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

2024-08-25 21:29:08: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2024-08-25 21:29:08: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

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: ['COPDEM30', 'TANDEM', 'DEM3', 'AW3D30', 'ARCTICDEM', 'ASTER', 'COPDEM90', 'ALASKA', 'MAPZEN']
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

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

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

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

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

df.corr()

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

df_diff.describe()

	COPDEM30-TANDEM	COPDEM30-DEM3	COPDEM30-AW3D30	COPDEM30-ARCTICDEM	COPDEM30-ASTER	COPDEM30-COPDEM90	COPDEM30-ALASKA	COPDEM30-MAPZEN	TANDEM-DEM3	TANDEM-AW3D30	...	ARCTICDEM-ASTER	ARCTICDEM-COPDEM90	ARCTICDEM-ALASKA	ARCTICDEM-MAPZEN	ASTER-COPDEM90	ASTER-ALASKA	ASTER-MAPZEN	COPDEM90-ALASKA	COPDEM90-MAPZEN	ALASKA-MAPZEN
count	1671.000000	1671.000000	0.0	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	1671.000000	0.0	...	767.000000	767.000000	767.000000	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	1671.000000	1671.000000
mean	123.708183	-51.034972	NaN	-110.319916	-18.073272	-0.180525	-39.221443	-41.711812	-174.743149	NaN	...	102.185113	110.184074	85.975136	83.946521	17.892747	-21.148171	-23.638540	-39.040920	-41.531287	-2.490369
std	144.721756	40.056016	NaN	114.278862	39.801988	5.457573	38.322636	38.776733	138.969170	NaN	...	108.948818	113.430321	113.401772	113.273900	38.488241	25.915689	25.367381	37.098961	37.537576	4.234587
min	-84.786133	-173.853027	NaN	-459.774902	-128.154297	-18.569824	-135.432617	-137.803711	-872.888672	NaN	...	-121.276367	-131.313965	-147.653809	-146.276367	-111.338867	-137.751953	-130.000000	-130.026367	-135.499023	-20.434082
25%	18.154541	-73.800049	NaN	-158.155029	-42.359863	-3.210938	-65.968750	-69.063965	-217.682861	NaN	...	28.157715	36.502441	5.880127	4.010010	-7.423828	-36.852539	-40.000000	-65.058350	-68.542236	-5.209473
50%	95.392578	-39.264160	NaN	-82.882324	-13.353027	-0.013672	-32.329102	-34.193848	-133.516602	NaN	...	71.489746	81.454102	58.945312	56.379883	12.738770	-19.398438	-22.000000	-33.313477	-34.726074	-2.204590
75%	178.104980	-22.958496	NaN	-35.619141	8.449463	2.751221	-10.385498	-13.209717	-84.589355	NaN	...	171.116211	158.221924	143.909668	143.172607	40.641357	-1.874023	-6.000000	-10.340332	-12.118164	-0.054199
max	812.099121	26.395508	NaN	136.359863	123.904297	17.572266	46.030762	38.260254	14.181152	NaN	...	428.646973	452.215332	401.538574	396.646973	124.032227	55.279785	54.000000	42.941895	33.757324	14.705078

8 rows × 36 columns

df_diff.abs().describe()

	COPDEM30-TANDEM	COPDEM30-DEM3	COPDEM30-AW3D30	COPDEM30-ARCTICDEM	COPDEM30-ASTER	COPDEM30-COPDEM90	COPDEM30-ALASKA	COPDEM30-MAPZEN	TANDEM-DEM3	TANDEM-AW3D30	...	ARCTICDEM-ASTER	ARCTICDEM-COPDEM90	ARCTICDEM-ALASKA	ARCTICDEM-MAPZEN	ASTER-COPDEM90	ASTER-ALASKA	ASTER-MAPZEN	COPDEM90-ALASKA	COPDEM90-MAPZEN	ALASKA-MAPZEN
count	1671.000000	1671.000000	0.0	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	1671.000000	0.0	...	767.000000	767.000000	767.000000	767.000000	1671.000000	1671.000000	1671.000000	1671.000000	1671.000000	1671.000000
mean	135.184494	52.023959	NaN	118.163139	32.801881	4.060298	42.954536	44.739207	174.815853	NaN	...	112.711556	117.711647	105.086456	103.447394	31.869296	26.224301	27.788151	42.339176	44.107472	3.799071
std	134.057739	38.762088	NaN	106.137970	28.886490	3.649907	34.083324	35.238571	138.877648	NaN	...	98.004051	105.587395	95.936180	95.769833	28.025463	20.760783	20.735849	33.283024	34.471963	3.113807
min	0.002441	0.012207	NaN	0.650879	0.005371	0.000488	0.007324	0.117188	0.848633	NaN	...	0.155762	0.624512	0.155273	0.289062	0.000000	0.001953	0.000000	0.000000	0.069336	0.002930
25%	38.805664	23.231689	NaN	41.187012	11.155762	1.359131	14.802979	16.361816	84.589355	NaN	...	41.791260	41.026611	36.557861	33.356445	10.204102	11.052246	12.000000	14.249756	14.985596	1.328125
50%	95.392578	39.264160	NaN	84.257812	23.466309	2.911621	32.889648	34.610840	133.516602	NaN	...	73.289062	82.970215	65.577148	64.257812	22.945801	21.200195	24.000000	33.496582	34.726074	3.072266
75%	178.104980	73.800049	NaN	158.155029	47.615723	5.748779	65.968750	69.063965	217.682861	NaN	...	171.116211	158.221924	144.205811	144.176514	46.310059	37.241211	40.000000	65.058350	68.542236	5.534912
max	812.099121	173.853027	NaN	459.774902	128.154297	18.569824	135.432617	137.803711	872.888672	NaN	...	428.646973	452.215332	401.538574	396.646973	124.032227	137.751953	130.000000	130.026367	135.499023	20.434082

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents