OGGM-Shop and Glacier Directories in OGGM

OGGM-Shop and Glacier Directories in OGGM#

New in OGGM 1.6: if you are using OGGM 1.6, all the datasets below are available as pre-processed directories. The tutorial below is useful if you want to add data from the shop to your existing directories.

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

2024-08-25 21:25:04: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:25:04: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:25:04: oggm.cfg: Multiprocessing: using all available processors (N=4)

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!

You can create a persistent OGGM working directory at a specific path via path = utils.mkdir(path). Beware! If you use reset=True in utils.mkdir, ALL DATA in this folder will be deleted!

Define the glaciers for the run#

rgi_ids = ['RGI60-01.13696']  # Malaspina glacier (large - hungry in memory)
rgi_ids = ['RGI60-05.00800']  #  (not all datasets)
rgi_ids = ['RGI60-14.06794']  # Baltoro glacier in Karakoram
rgi_ids = ['RGI60-11.00897']  # Hintereisferner (not all datasets)
rgi_ids = ['RGI60-01.10689']  # Columbia glacier

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

2024-08-25 21:25:06: oggm.workflow: init_glacier_directories from prepro level 1 on 1 glaciers.
2024-08-25 21:25:06: oggm.workflow: Execute entity tasks [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]
print('Path to the DEM:', gdir.get_filepath('glacier_mask'))

Path to the DEM: /tmp/OGGM/OGGM-Shop/per_glacier/RGI60-01/RGI60-01.10/RGI60-01.10689/glacier_mask.tif

And we can also access some attributes of this glacier:

gdir

<oggm.GlacierDirectory>
  RGI id: RGI60-01.10689
  Region: 01: Alaska
  Subregion: 01-04: W Chugach Mtns (Talkeetna)      
  Name: Columbia Glacier
  Glacier type: Glacier
  Terminus type: Marine-terminating
  Status: Glacier or ice cap
  Area: 773.873 km2
  Lon, Lat: (-146.888, 61.299)
  Grid (nx, ny): (223, 233)
  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))

../../_images/62a780504538de3cf1b598f7375c21ccb821390640c6fdc4c267aa36876a21d2.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-01.10689
Available DEM sources: ['COPDEM30', 'TANDEM', 'DEM3', 'AW3D30', 'ARCTICDEM', 'ASTER', 'COPDEM90', 'ALASKA', 'MAPZEN']

# We use xarray to store the data
import rioxarray as rioxr
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 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)
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/d1dfbbe14391283d3e116a25fd6bec6d556d18dbd4dec79e9853c464ee32a05f.png

Original (raw) topography data#

See dem_sources.ipynb.

OGGM-Shop: bed topography data from Farinotti et al., (2019)#

OGGM can also download data from the Farinotti et al., (2019) consensus estimate and reproject it to the glacier directories map:

from oggm.shop import bedtopo

workflow.execute_entity_task(bedtopo.add_consensus_thickness, gdirs);

2024-08-25 21:25:17: oggm.workflow: Execute entity tasks [add_consensus_thickness] on 1 glaciers

with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    ds = ds.load()

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)');
ax.set_title('Farinotti 19 thickness');

../../_images/a0a1ca53a3b0e87ffa37dd5b7920e4859d7211e0e2586fd0721d94d3a77025f2.png

OGGM-Shop: climate data#

# TODO

What’s next?#

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