OGGM-Shop and Glacier Directories in OGGM

Set-up

Input data folders

If you are using your own computer: before you start, make sure that you have set-up the input data configuration file at your wish. In the course of this tutorial, we will need to download data needed for each glacier (a couple of mb at max, depending on the chosen glaciers), so make sure you have an internet connection.

cfg.initialize() and cfg.PARAMS

An OGGM simulation script will always start with the following commands:

from oggm import cfg, utils
cfg.initialize(logging_level='WARNING')
2021-05-18 08:42:36: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2021-05-18 08:42:36: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2021-05-18 08:42:36: oggm.cfg: Multiprocessing: using all available processors (N=2)

A call to cfg.initialize() will read the default parameter file (or any user-provided file) and make them available to all other OGGM tools via the cfg.PARAMS dictionary. Here are some examples of these parameters:

cfg.PARAMS['continue_on_error']
False
cfg.PARAMS['border']
40
cfg.PARAMS['has_internet']
True

Workflow

import os
from oggm import workflow, tasks

Working directory

Each OGGM run needs a single folder where to store the results of the computations for all glaciers. This is called a “working directory” and needs to be specified before each run. Here we create a temporary folder for you:

cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-Shop', reset=True)
cfg.PATHS['working_dir']
'/tmp/OGGM/OGGM-Shop'

We use a temporary directory for this example, but in practice you will set this working directory yourself (for example: /home/john/OGGM_output. The size of this directory will depend on how many glaciers you’ll simulate!

This working directory is meant to be persistent, i.e. you can stop your processing workflow after any task, and restart from an existing working directory at a later stage.

Define the glaciers for the run

rgi_ids = ['RGI60-01.13696']  # Malaspina glacier (large - hungry in memory)
rgi_ids = ['RGI60-05.00800']  # Glacier in Greenland

You can provide any number of glacier identifiers. You can find other glacier identifiers by exploring the GLIMS viewer.

For an operational run on an RGI region, you might want to download the Randolph Glacier Inventory dataset instead, and start from it. This case is covered in the working with the RGI tutorial and in the “Starting from scratch” section below.

Starting from RGItopo

The OGGM workflow is organized as a list of tasks that have to be applied to a list of glaciers. The vast majority of tasks are called entity tasks: they are standalone operations to be realized on one single glacier entity. These tasks are executed sequentially (one after another): they often need input generated by the previous task(s): for example, the glacier mask task needs the glacier topography data.

To handle this situation, OGGM uses a workflow based on data persistence on disk: instead of passing data as python variables from one task to another, each task will read the data from disk and then write the computation results back to the disk, making these new data available for the next task in the queue.

These glacier specific data are located in glacier directories. These directories are initialized with the following command (this can take a little while on the first call, as OGGM needs to download some data):

# 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.
prepro_border = 10
# Degree of processing level. This is OGGM specific and for the shop 1 is the one you want
from_prepro_level = 1
# URL of the preprocessed Gdirs
base_url = 'https://cluster.klima.uni-bremen.de/data/gdirs/dems_v1/default/'

gdirs = workflow.init_glacier_directories(rgi_ids,
                                          from_prepro_level=from_prepro_level,
                                          prepro_base_url=base_url,
                                          prepro_rgi_version=prepro_rgi_version,
                                          prepro_border=prepro_border)
2021-05-18 08:42:37: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2021-05-18 08:42:37: oggm.workflow: Execute entity task gdir_from_prepro on 1 glaciers

gdirs is a list of GlacierDirectory objects (one for each glacier). Glacier directories are used by OGGM as “file and attribute manager” for single glaciers.

For example, we now know where to find the glacier mask files for this glacier:

gdir = gdirs[0]  # take Unteraar
print('Path to the DEM:', gdir.get_filepath('glacier_mask'))
Path to the DEM: /tmp/OGGM/OGGM-Shop/per_glacier/RGI60-05/RGI60-05.00/RGI60-05.00800/glacier_mask.tif

And we can also access some attributes of this glacier:

gdir
<oggm.GlacierDirectory>
  RGI id: RGI60-05.00800
  Region: 05: Greenland
  Subregion: 05-01: Greenland (periphery)           
  Glacier type: Glacier
  Terminus type: Marine-terminating
  Status: Glacier or ice cap
  Area: 246.825 km2
  Lon, Lat: (-52.141, 66.0945)
  Grid (nx, ny): (176, 158)
  Grid (dx, dy): (200.0, -200.0)
gdir.rgi_date  # date at which the outlines are valid
2000

The advantage of this Glacier Directory data model is that it simplifies greatly the data transfer between tasks. The single mandatory argument of all entity tasks will allways be a glacier directory. With the glacier directory, each task will find the input it needs: for example, the glacier outlines are needed for the next plotting function, and are available via the gdir argument:

from oggm import graphics
graphics.plot_googlemap(gdir, figsize=(8, 7))
../_images/oggm_shop_29_0.png

For most glaciers in the world there are several digital elevation models (DEM) which cover the respective glacier. In OGGM we have currently implemented many different open access DEMs to choose from. For some, you need to register to get access, see dem_sources.ipynb/register. 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.

RGItopo data

sources = [src for src in os.listdir(gdir.dir) if src in utils.DEM_SOURCES]
print('RGI ID:', gdir.rgi_id)
print('Available DEM sources:', sources)
RGI ID: RGI60-05.00800
Available DEM sources: ['ASTER', 'AW3D30', 'COPDEM', 'ARCTICDEM', 'MAPZEN', 'DEM3', 'TANDEM', 'GIMP']
# We use xarray to store the data
import xarray as xr
import numpy as np

ods = xr.Dataset()
for src in sources:
    demfile = os.path.join(gdir.dir, src) + '/dem.tif'
    with xr.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 xr.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)
n_col = 3
x_size = 12
n_rows = -(-ns // n_col)
y_size = x_size / n_col * n_rows
from mpl_toolkits.axes_grid1 import AxesGrid
import salem
import matplotlib.pyplot as plt

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_col),
                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]:
    grid[-1].remove()
    grid[-1].cax.remove()
../_images/oggm_shop_36_0.png

Original (raw) topography data

See dem_sources.ipynb.

OGGM-Shop: ITS-live

This is an example on how to extract velocity fields from the ITS_live Regional Glacier and Ice Sheet Surface Velocities Mosaic (Gardner, A. et al 2019) at 120 m resolution and reproject this data to the OGGM-glacier grid. This only works where ITS-live data is available! (not in the Alps)

# this will download severals large dataset (2*~800MB)
from oggm.shop import its_live, rgitopo
workflow.execute_entity_task(rgitopo.select_dem_from_dir, gdirs, dem_source='COPDEM', keep_dem_folders=True);
workflow.execute_entity_task(tasks.glacier_masks, gdirs);
workflow.execute_entity_task(its_live.velocity_to_gdir, gdirs);
2021-05-18 08:43:01: oggm.workflow: Execute entity task select_dem_from_dir on 1 glaciers
2021-05-18 08:43:01: oggm.workflow: Execute entity task glacier_masks on 1 glaciers
2021-05-18 08:43:03: oggm.workflow: Execute entity task velocity_to_gdir on 1 glaciers

By applying the entity task its_live.velocity_to_gdir() the model downloads and reprojects the ITS_live files to a given glacier map.

The velocity components (vx, vy) are added to the gridded_data nc file stored on each glacier directory.

According to the ITS_LIVE documentation velocities are given in ground units (i.e. absolute velocities). We then use bilinear interpolation to reproject the velocities to the local glacier map by re-projecting the vector distances.

By specifying add_error=True, we also reproject and scale the error for each component (evx, evy).

Now we can read in all the gridded data that comes with OGGM, including the ITS_Live velocity components.

with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    ds = ds.load()
ds
<xarray.Dataset>
Dimensions:          (x: 176, y: 158)
Coordinates:
  * x                (x) float32 -1.828e+04 -1.808e+04 ... 1.652e+04 1.672e+04
  * y                (y) float32 7.344e+06 7.344e+06 ... 7.313e+06 7.312e+06
Data variables:
    topo             (y, x) float32 1.56e+03 1.576e+03 1.596e+03 ... 675.4 634.6
    topo_smoothed    (y, x) float32 1.562e+03 1.576e+03 ... 666.6 640.1
    topo_valid_mask  (y, x) int8 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1
    glacier_mask     (y, x) int8 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0
    glacier_ext      (y, x) int8 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0
    obs_icevel_x     (y, x) float32 0.1185 0.1567 0.1677 ... 7.448 5.996 4.791
    obs_icevel_y     (y, x) float32 1.509 0.8505 0.5083 ... 7.928 5.952 3.744
Attributes:
    author:         OGGM
    author_info:    Open Global Glacier Model
    pyproj_srs:     +proj=tmerc +lat_0=0 +lon_0=-52.141 +k=0.9996 +x_0=0 +y_0...
    max_h_dem:      2159.736
    min_h_dem:      0.0
    max_h_glacier:  2097.7734
    min_h_glacier:  0.0
# plot the salem map background, make countries in grey
smap = ds.salem.get_map(countries=False)
smap.set_shapefile(gdir.read_shapefile('outlines'))
smap.set_topography(ds.topo.data);
# get the velocity data
u = ds.obs_icevel_x.where(ds.glacier_mask == 1)
v = ds.obs_icevel_y.where(ds.glacier_mask == 1)
ws = (u**2 + v**2)**0.5

The ds.glacier_mask == 1 command will remove the data outside of the glacier outline.

# get the axes ready
f, ax = plt.subplots(figsize=(9, 9))

# Quiver only every 3rd grid point
us = u[1::3, 1::3]
vs = v[1::3, 1::3]

smap.set_data(ws)
smap.set_cmap('Blues')
smap.plot(ax=ax)
smap.append_colorbar(ax=ax, label = 'ice velocity (m yr$^{-1}$)')

# transform their coordinates to the map reference system and plot the arrows
xx, yy = smap.grid.transform(us.x.values, us.y.values, crs=gdir.grid.proj)
xx, yy = np.meshgrid(xx, yy)
qu = ax.quiver(xx, yy, us.values, vs.values)
qk = ax.quiverkey(qu, 0.82, 0.97, 100, '100 m yr$^{-1}$',
                   labelpos='E', coordinates='axes')
../_images/oggm_shop_47_0.png

OGGM-Shop: bed topography data

from oggm.shop import bedtopo
workflow.execute_entity_task(bedtopo.add_consensus_thickness, gdirs);
2021-05-18 08:44:07: oggm.workflow: Execute entity task add_consensus_thickness on 1 glaciers
with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    ds = ds.load()
# plot the salem map background, make countries in grey
smap = ds.salem.get_map(countries=False)
smap.set_shapefile(gdir.read_shapefile('outlines'))
smap.set_topography(ds.topo.data);
f, ax = plt.subplots(figsize=(9, 9))

smap.set_data(ds.consensus_ice_thickness)
smap.set_cmap('Blues')
smap.plot(ax=ax)
smap.append_colorbar(ax=ax, label='ice thickness (m)');
../_images/oggm_shop_53_0.png

OGGM-Shop: climate data

# TODO

What’s next?