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

2023-08-28 00:06:22: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-08-28 00:06:22: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-08-28 00:06:22: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-08-28 00:06:23: 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]

2023-08-28 00:06:24: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2023-08-28 00:06:24: 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: Land-terminating
  Status: Glacier or ice cap
  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: ['AW3D30', 'COPDEM90', 'DEM3', 'ASTER', 'MAPZEN', 'COPDEM30', 'ARCTICDEM', 'ALASKA', 'TANDEM']
Plotting directory: /__w/tutorials/tutorials/notebooks/others/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

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

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

/usr/local/pyenv/versions/3.10.12/lib/python3.10/site-packages/seaborn/axisgrid.py:118: UserWarning: The figure layout has changed to tight
  self._figure.tight_layout(*args, **kwargs)

Table statistics

df.describe()

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

df.corr()

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

df_diff.describe()

	AW3D30-COPDEM90	AW3D30-DEM3	AW3D30-ASTER	AW3D30-MAPZEN	AW3D30-COPDEM30	AW3D30-ARCTICDEM	AW3D30-ALASKA	AW3D30-TANDEM	COPDEM90-DEM3	COPDEM90-ASTER	...	MAPZEN-COPDEM30	MAPZEN-ARCTICDEM	MAPZEN-ALASKA	MAPZEN-TANDEM	COPDEM30-ARCTICDEM	COPDEM30-ALASKA	COPDEM30-TANDEM	ARCTICDEM-ALASKA	ARCTICDEM-TANDEM	ALASKA-TANDEM
count	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1671.000000	1671.000000	...	1671.000000	767.000000	1671.000000	1671.000000	767.000000	1671.000000	1671.000000	767.000000	767.000000	1671.000000
mean	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-50.854447	-17.892747	...	41.711812	-83.946521	2.490369	165.419989	-110.319916	-39.221443	123.708183	85.975136	251.653717	162.929626
std	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	38.750864	38.488241	...	38.776733	113.273900	4.234587	142.198023	114.278862	38.322636	144.721756	113.401764	224.423676	142.417313
min	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-163.009766	-124.032227	...	-38.260254	-396.646973	-14.705078	-16.507324	-459.774902	-135.432617	-84.786133	-147.653809	-48.799805	-15.732910
25%	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-74.026611	-40.641357	...	13.209717	-143.172607	0.054199	71.849121	-158.155029	-65.968750	18.154541	5.880127	111.032227	67.788818
50%	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-39.517090	-12.738770	...	34.193848	-56.379883	2.204590	123.718750	-82.882324	-32.329102	95.392578	58.945312	170.058594	120.674805
75%	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	-23.536377	7.423828	...	69.063965	-4.010010	5.209473	214.740479	-35.619141	-10.385498	178.104980	143.909668	358.547607	212.847412
max	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	17.004395	111.338867	...	137.803711	146.276367	20.434082	851.888672	136.359863	46.030762	812.099121	401.538574	1085.924805	848.312500

8 rows × 36 columns

df_diff.abs().describe()

	AW3D30-COPDEM90	AW3D30-DEM3	AW3D30-ASTER	AW3D30-MAPZEN	AW3D30-COPDEM30	AW3D30-ARCTICDEM	AW3D30-ALASKA	AW3D30-TANDEM	COPDEM90-DEM3	COPDEM90-ASTER	...	MAPZEN-COPDEM30	MAPZEN-ARCTICDEM	MAPZEN-ALASKA	MAPZEN-TANDEM	COPDEM30-ARCTICDEM	COPDEM30-ALASKA	COPDEM30-TANDEM	ARCTICDEM-ALASKA	ARCTICDEM-TANDEM	ALASKA-TANDEM
count	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1671.000000	1671.000000	...	1671.000000	767.000000	1671.000000	1671.000000	767.000000	1671.000000	1671.000000	767.000000	767.000000	1671.000000
mean	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	51.431971	31.869296	...	44.739207	103.447394	3.799071	165.650691	118.163139	42.954536	135.184494	105.086456	253.502304	163.309280
std	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	37.980538	28.025463	...	35.238571	95.769833	3.113807	141.929044	106.137970	34.083321	134.057739	95.936180	222.330795	141.981552
min	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.002441	0.000000	...	0.117188	0.289062	0.002930	0.151367	0.650879	0.007324	0.002441	0.155273	0.293457	0.070801
25%	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	23.536377	10.204102	...	16.361816	33.356445	1.328125	71.849121	41.187012	14.802979	38.805664	36.557861	111.032227	67.788818
50%	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	39.517090	22.945801	...	34.610840	64.257812	3.072266	123.718750	84.257812	32.889648	95.392578	65.577148	170.058594	120.674805
75%	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	74.026611	46.310059	...	69.063965	144.176514	5.534912	214.740479	158.155029	65.968750	178.104980	144.205811	358.547607	212.847412
max	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	163.009766	124.032227	...	137.803711	396.646973	20.434082	851.888672	459.774902	135.432617	812.099121	401.538574	1085.924805	848.312500

8 rows × 36 columns

What’s next?

• return to the OGGM documentation

• back to the table of contents