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-02-17 22:20:46: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.

2021-02-17 22:20:46: oggm.cfg: Multiprocessing switched OFF according to the parameter file.

2021-02-17 22:20:46: oggm.cfg: Multiprocessing: using all available processors (N=8)

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

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-02-17 22:20:47: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.

2021-02-17 22:20:47: 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
  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

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

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: ['DEM3', 'ASTER', 'GIMP', 'COPDEM', 'ARCTICDEM', 'TANDEM', 'MAPZEN', 'AW3D30']

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

Original (raw) topography data¶

See dem_sources.ipynb.

OGGM-Shop: bed topography data¶

from oggm.shop import bedtopo

workflow.execute_entity_task(bedtopo.add_consensus_thickness, gdirs);

2021-02-17 22:21:18: 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)');

OGGM-Shop: climate data¶

# TODO

What’s next?¶

look at the OGGM-Shop documentation
return to the OGGM documentation
back to the table of contents

OGGM tutorials