1. Preprocess a subset of an RGI Region

This example shows how to run the first steps of the OGGM preprocessing chain for a subset of the Alps (the Rofental catchment in the Austrian Alps).

Script

# Python imports
import os

# Libs
import geopandas as gpd
import shapely.geometry as shpg

# Locals
import oggm.cfg as cfg
from oggm import utils, workflow

# For timing the run
import time
start = time.time()

# Initialize OGGM and set up the default run parameters
cfg.initialize()
rgi_version = '61'
rgi_region = '11'  # Alps

# Local working directory (where OGGM will write its output)
WORKING_DIR = utils.gettempdir('OGGM_Rofental')
utils.mkdir(WORKING_DIR, reset=True)
cfg.PATHS['working_dir'] = WORKING_DIR

# We use intersects
path = utils.get_rgi_intersects_region_file(rgi_region, version=rgi_version)
cfg.set_intersects_db(path)

# RGI file
path = utils.get_rgi_region_file(rgi_region, version=rgi_version)
rgidf = gpd.read_file(path)

# Get the Rofental Basin file
path = utils.get_demo_file('rofental_hydrosheds.shp')
basin = gpd.read_file(path)

# Take all glaciers in the Rofental Basin
in_bas = [basin.geometry.contains(shpg.Point(x, y))[0] for
          (x, y) in zip(rgidf.CenLon, rgidf.CenLat)]
rgidf = rgidf.loc[in_bas]
# Store them for later
rgidf.to_file(os.path.join(WORKING_DIR, 'rgi_rofental.shp'))

# Sort for more efficient parallel computing
rgidf = rgidf.sort_values('Area', ascending=False)

print('Starting OGGM run')
print('Number of glaciers: {}'.format(len(rgidf)))

# Go - initialize glacier directories
gdirs = workflow.init_glacier_regions(rgidf)

# Tasks shortcuts - see the next examples for more details
workflow.gis_prepro_tasks(gdirs)
workflow.climate_tasks(gdirs)
workflow.inversion_tasks(gdirs)

# Compile output
print('Compiling output')
utils.compile_glacier_statistics(gdirs)
utils.write_centerlines_to_shape(gdirs)

# Log
m, s = divmod(time.time() - start, 60)
h, m = divmod(m, 60)
print('OGGM is done! Time needed: %d:%02d:%02d' % (h, m, s))

If everything went well, you should see an output similar to:

2018-01-19 16:56:03: oggm.cfg: Parameter file:
Starting OGGM run
Number of glaciers: 54
2018-01-19 16:56:04: oggm.workflow: Multiprocessing: using all available processors (N=8)
2018-01-19 16:56:04: oggm.core.gis: (RGI60-11.00719) define_glacier_region
(...)
2018-01-19 16:56:20: oggm.core.inversion: (RGI60-11.00873) filter_inversion_output
2018-01-19 16:56:20: oggm.core.inversion: (RGI60-11.00747) filter_inversion_output
2018-01-19 16:56:20: oggm.core.inversion: (RGI60-11.01040) filter_inversion_output
Compiling output
OGGM is done! Time needed: 0:00:19

Some analyses

The output directory contains some compiled output files. The compile_glacier_statistics.csv file contains various information for each glacier after the preprocessing, and we also wrote the centerlines as a shapefile. Let’s plot it:

# Imports
from os import path
import geopandas as gpd
import matplotlib.pyplot as plt
from oggm.utils import get_demo_file, gettempdir

# Local working directory (where OGGM wrote its output)
WORKING_DIR = gettempdir('OGGM_Rofental')

# Plot: the basin, the outlines and the centerlines
basin = gpd.read_file(get_demo_file('rofental_hydrosheds.shp'))
rgi = gpd.read_file(path.join(WORKING_DIR, 'rgi_rofental.shp'))
centerlines = gpd.read_file(path.join(WORKING_DIR, 'glacier_centerlines.shp'))

f, ax = plt.subplots()
basin.plot(ax=ax, color='k', alpha=0.2)
rgi.plot(ax=ax, color='C0')
centerlines.plot(ax=ax, color='C3')
plt.title('Rofental glaciers and centerlines')
plt.tight_layout()
plt.show()

This code snippet should produce the following plot:

../_images/rofental_centerlines.png