Ice thickness inversion
Contents
Ice thickness inversion#
This example shows how to run the OGGM ice thickness inversion model with various ice parameters: the deformation parameter A and a sliding parameter (fs).
There is currently no “best” set of parameters for the ice thickness inversion model. As shown in Maussion et al. (2019), the default parameter set results in global volume estimates which are a bit larger than previous values. For the consensus estimate of Farinotti et al. (2019), OGGM participated with a deformation parameter A that is 1.5 times larger than the generally accepted default value.
There is no reason to think that the ice parameters are the same between neighboring glaciers. There is currently no “good” way to calibrate them, or at least no generaly accepted one. We won’t discuss the details here, but we provide a script to illustrate the sensitivity of the model to this choice.
New in version 1.4: we demonstrate how to apply a new global task in OGGM, workflow.calibrate_inversion_from_consensus to calibrate the A parameter to match the consensus estimate from Farinotti et al. (2019).
At the end of this tutorial, we show how to distribute the “flowline thicknesses” on a glacier map.
Run#
# Libs
import geopandas as gpd
# Locals
import oggm.cfg as cfg
from oggm import utils, workflow, tasks, graphics
# Initialize OGGM and set up the default run parameters
cfg.initialize(logging_level='WARNING')
rgi_region = '11' # Region Central Europe
# Local working directory (where OGGM will write its output)
WORKING_DIR = utils.gettempdir('OGGM_Inversion')
cfg.PATHS['working_dir'] = WORKING_DIR
# This is useful here
cfg.PARAMS['use_multiprocessing'] = True
# RGI file
path = utils.get_rgi_region_file(rgi_region)
rgidf = gpd.read_file(path)
# Select the glaciers in the Pyrenees
rgidf = rgidf.loc[rgidf['O2Region'] == '2']
# Sort for more efficient parallel computing
rgidf = rgidf.sort_values('Area', ascending=False)
# Go - get the pre-processed glacier directories
# We start at level 3, because we need all data for the inversion
gdirs = workflow.init_glacier_directories(rgidf, from_prepro_level=3, prepro_border=10)
# Because of recent changes in OGGM not yet available in the preprocessed directories,
# we re-run this task:
workflow.execute_entity_task(tasks.prepare_for_inversion, gdirs)
# Default parameters
# Deformation: from Cuffey and Patterson 2010
glen_a = 2.4e-24
# Sliding: from Oerlemans 1997
fs = 5.7e-20
2023-03-07 12:29:35: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-07 12:29:35: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-07 12:29:35: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-03-07 12:29:35: oggm.cfg: Multiprocessing switched ON after user settings.
2023-03-07 12:29:46: oggm.workflow: init_glacier_directories from prepro level 3 on 35 glaciers.
2023-03-07 12:29:46: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 35 glaciers
2023-03-07 12:29:51: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
with utils.DisableLogger(): # this scraps some output - to use with caution!!!
# Correction factors
factors = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
factors += [1.1, 1.2, 1.3, 1.5, 1.7, 2, 2.5, 3, 4, 5]
factors += [6, 7, 8, 9, 10]
# Run the inversions tasks with the given factors
for f in factors:
# Without sliding
suf = '_{:03d}_without_fs'.format(int(f * 10))
workflow.execute_entity_task(tasks.mass_conservation_inversion, gdirs,
glen_a=glen_a*f, fs=0)
# Store the results of the inversion only
utils.compile_glacier_statistics(gdirs, filesuffix=suf,
inversion_only=True)
# With sliding
suf = '_{:03d}_with_fs'.format(int(f * 10))
workflow.execute_entity_task(tasks.mass_conservation_inversion, gdirs,
glen_a=glen_a*f, fs=fs)
# Store the results of the inversion only
utils.compile_glacier_statistics(gdirs, filesuffix=suf,
inversion_only=True)
Read the data#
The data are stored as csv files in the working directory. The easiest way to read them is to use pandas!
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy import stats
import os
Let’s read the output of the inversion with the default OGGM parameters first:
df = pd.read_csv(os.path.join(WORKING_DIR, 'glacier_statistics_011_without_fs.csv'), index_col=0)
df
| rgi_region | rgi_subregion | name | cenlon | cenlat | rgi_area_km2 | rgi_year | glacier_type | terminus_type | is_tidewater | ... | inv_volume_km3 | vas_volume_km3 | inv_volume_bsl_km3 | inv_volume_bwl_km3 | dem_source | n_orig_centerlines | flowline_type | perc_invalid_flowline | inversion_glen_a | inversion_fs | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| rgi_id | |||||||||||||||||||||
| RGI60-11.03208 | 11 | 11-02 | Aneto | 0.646032 | 42.641357 | 0.622 | 2011 | Glacier | Land-terminating | False | ... | 0.029854 | 0.017699 | 0.0 | 0.0 | NASADEM | 2 | centerlines | 0.205128 | 2.640000e-24 | 0 |
| RGI60-11.03232 | 11 | 11-02 | Ossoue | -0.141222 | 42.771191 | 0.449 | 2011 | Glacier | Land-terminating | False | ... | 0.019636 | 0.011306 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03222 | 11 | 11-02 | Monte Perdido Lower | 0.039655 | 42.679337 | 0.308 | 2011 | Glacier | Land-terminating | False | ... | 0.009120 | 0.006733 | 0.0 | 0.0 | NASADEM | 2 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03209 | 11 | 11-02 | Maladeta E | 0.639343 | 42.648979 | 0.260 | 2011 | Glacier | Land-terminating | False | ... | 0.007672 | 0.005334 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03231 | 11 | 11-02 | Oulettes Du Gaube | -0.144215 | 42.779716 | 0.130 | 2011 | Glacier | Land-terminating | False | ... | 0.002747 | 0.002057 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03213 | 11 | 11-02 | Seil De La Baque E | 0.494304 | 42.694481 | 0.105 | 2011 | Glacier | Land-terminating | False | ... | 0.002118 | 0.001533 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03239 | 11 | 11-02 | Jezerces Gl No1 | 19.807627 | 42.465839 | 0.097 | 2007 | Glacier | Land-terminating | False | ... | 0.002538 | 0.001375 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.562500 | 2.640000e-24 | 0 |
| RGI60-11.03229 | 11 | 11-02 | Gabietous | -0.057456 | 42.695129 | 0.087 | 2011 | Glacier | Land-terminating | False | ... | 0.001720 | 0.001184 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03218 | 11 | 11-02 | Llardana | 0.428193 | 42.653946 | 0.084 | 2011 | Glacier | Land-terminating | False | ... | 0.001472 | 0.001128 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03228 | 11 | 11-02 | Taillon | -0.039683 | 42.695419 | 0.083 | 2011 | Glacier | Land-terminating | False | ... | 0.001733 | 0.001110 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03205 | 11 | 11-02 | Tempestades | 0.666092 | 42.626839 | 0.077 | 2011 | Glacier | Land-terminating | False | ... | 0.001395 | 0.001001 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03206 | 11 | 11-02 | Barrancs | 0.658785 | 42.632011 | 0.067 | 2011 | Glacier | Land-terminating | False | ... | 0.001078 | 0.000827 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03217 | 11 | 11-02 | La Paul | 0.437356 | 42.661072 | 0.057 | 2011 | Glacier | Land-terminating | False | ... | 0.000993 | 0.000662 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03221 | 11 | 11-02 | Monte Perdido Upper | 0.033972 | 42.678589 | 0.054 | 2011 | Glacier | Land-terminating | False | ... | 0.001139 | 0.000614 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03210 | 11 | 11-02 | Maladeta W | 0.632566 | 42.652107 | 0.054 | 2011 | Glacier | Land-terminating | False | ... | 0.000884 | 0.000614 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03233 | 11 | 11-02 | Infierno Central | -0.258849 | 42.783802 | 0.053 | 2011 | Glacier | Land-terminating | False | ... | 0.000868 | 0.000599 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03225 | 11 | 11-02 | Astazou | 0.034118 | 42.703136 | 0.052 | 2011 | Glacier | Land-terminating | False | ... | 0.001376 | 0.000583 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03211 | 11 | 11-02 | Boum | 0.555278 | 42.700157 | 0.050 | 2011 | Glacier | Land-terminating | False | ... | 0.001101 | 0.000553 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03240 | 11 | 11-02 | Jezerces Gl No2 | 19.815450 | 42.447353 | 0.049 | 2007 | Glacier | Land-terminating | False | ... | 0.000666 | 0.000538 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03227 | 11 | 11-02 | Pailla W | 0.015287 | 42.705482 | 0.046 | 2011 | Glacier | Land-terminating | False | ... | 0.000711 | 0.000493 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03219 | 11 | 11-02 | Barroude | 0.142224 | 42.726147 | 0.042 | 2011 | Glacier | Land-terminating | False | ... | 0.000828 | 0.000435 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03220 | 11 | 11-02 | Munia | 0.129816 | 42.717735 | 0.041 | 2011 | Glacier | Land-terminating | False | ... | 0.000583 | 0.000421 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03212 | 11 | 11-02 | Portillon D Oo | 0.511460 | 42.693188 | 0.040 | 2011 | Glacier | Land-terminating | False | ... | 0.000549 | 0.000407 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03226 | 11 | 11-02 | Pailla E | 0.021260 | 42.705322 | 0.035 | 2011 | Glacier | Land-terminating | False | ... | 0.000639 | 0.000339 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03241 | 11 | 11-02 | It4L30000001 Gh Del Calderone | 13.566710 | 42.470540 | 0.033 | 1990 | Glacier | Land-terminating | False | ... | 0.000619 | 0.000312 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.250000 | 2.640000e-24 | 0 |
| RGI60-11.03230 | 11 | 11-02 | Petit Vignemale | -0.137718 | 42.775856 | 0.027 | 2011 | Glacier | Land-terminating | False | ... | 0.000370 | 0.000237 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03215 | 11 | 11-02 | Gourgs Blancs | 0.480293 | 42.700859 | 0.026 | 2011 | Glacier | Land-terminating | False | ... | 0.000468 | 0.000225 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03234 | 11 | 11-02 | Las Neous | -0.286945 | 42.838642 | 0.025 | 2011 | Glacier | Land-terminating | False | ... | 0.000465 | 0.000213 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03223 | 11 | 11-02 | Marbore | 0.018435 | 42.693417 | 0.023 | 2011 | Glacier | Land-terminating | False | ... | 0.000325 | 0.000190 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03224 | 11 | 11-02 | Cilindro | 0.026603 | 42.689346 | 0.022 | 2011 | Glacier | Land-terminating | False | ... | 0.000319 | 0.000179 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03204 | 11 | 11-02 | Valier | 1.088501 | 42.798626 | 0.021 | 2011 | Glacier | Land-terminating | False | ... | 0.000264 | 0.000168 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03214 | 11 | 11-02 | Seil De La Baque W | 0.490498 | 42.696220 | 0.020 | 2011 | Glacier | Land-terminating | False | ... | 0.000303 | 0.000157 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03238 | 11 | 11-02 | Debeli Namet Gl | 19.067500 | 43.115000 | 0.019 | 2003 | Glacier | Land-terminating | False | ... | 0.000180 | 0.000146 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03207 | 11 | 11-02 | Coronas | 0.654112 | 42.632977 | 0.015 | 2011 | Glacier | Land-terminating | False | ... | 0.000259 | 0.000106 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
| RGI60-11.03216 | 11 | 11-02 | Posets | 0.438703 | 42.656570 | 0.010 | 2011 | Glacier | Land-terminating | False | ... | 0.000130 | 0.000060 | 0.0 | 0.0 | NASADEM | 1 | centerlines | 0.000000 | 2.640000e-24 | 0 |
35 rows × 21 columns
There are only 35 glaciers in the Pyrenees! That’s why the run was relatively fast.
Results#
One way to visualize the output is to plot the volume as a function of area in a log-log plot, illustrating the well known volume-area relationship of mountain glaciers:
ax = df.plot(kind='scatter', x='rgi_area_km2', y='inv_volume_km3')
ax.semilogx(); ax.semilogy()
xlim, ylim = [1e-2, 0.7], [1e-5, 0.05]
ax.set_xlim(xlim); ax.set_ylim(ylim);
As we can see, there is a clear relationship, but it is not perfect. Let’s fit a line to these data (the “volume-area scaling law”):
# Fit in log space
dfl = np.log(df[['inv_volume_km3', 'rgi_area_km2']])
slope, intercept, r_value, p_value, std_err = stats.linregress(dfl.rgi_area_km2.values, dfl.inv_volume_km3.values)
In their seminal paper, Bahr et al. (1997) describe this relationship as:
\[V = \alpha S^{\gamma}\]
With V the volume in km\(^3\), S the area in km\(^2\) and \(\alpha\) and \(\gamma\) the scaling parameters (0.034 and 1.375, respectively). How does OGGM compare to these in the Pyrenees?
print('power: {:.3f}'.format(slope))
print('slope: {:.3f}'.format(np.exp(intercept)))
power: 1.301
slope: 0.045
ax = df.plot(kind='scatter', x='rgi_area_km2', y='inv_volume_km3', label='OGGM glaciers')
ax.plot(xlim, np.exp(intercept) * (xlim ** slope), color='C3', label='Fitted line')
ax.semilogx(); ax.semilogy()
ax.set_xlim(xlim); ax.set_ylim(ylim);
ax.legend();
Sensitivity analysis#
Now, let’s read the output files of each run separately, and compute the regional volume out of them:
dftot = pd.DataFrame(index=factors)
for f in factors:
# Without sliding
suf = '_{:03d}_without_fs'.format(int(f * 10))
fpath = os.path.join(WORKING_DIR, 'glacier_statistics{}.csv'.format(suf))
_df = pd.read_csv(fpath, index_col=0, low_memory=False)
dftot.loc[f, 'without_sliding'] = _df.inv_volume_km3.sum()
# With sliding
suf = '_{:03d}_with_fs'.format(int(f * 10))
fpath = os.path.join(WORKING_DIR, 'glacier_statistics{}.csv'.format(suf))
_df = pd.read_csv(fpath, index_col=0, low_memory=False)
dftot.loc[f, 'with_sliding'] = _df.inv_volume_km3.sum()
And plot them:
dftot.plot();
plt.xlabel('Factor of Glen A (default 1)'); plt.ylabel('Regional volume (km$^3$)');
As you can see, there is quite a difference between the solutions. In particular, close to the default value for Glen A, the regional estimates are very sensitive to small changes in A. The calibration of A is a problem that has yet to be resolved by global glacier models…
New in version 1.4: calibrate to match the consensus estimate#
Here, a “best Glen A” is found in order that the total inverted volume of the glaciers of gdirs fits to the 2019 consensus estimate.
cdf = workflow.calibrate_inversion_from_consensus(gdirs, filter_inversion_output=False)
2023-03-07 12:35:07: oggm.workflow: Applying global task calibrate_inversion_from_consensus on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Consensus estimate optimisation with A factor: 0.1 and fs: 0
2023-03-07 12:35:08: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [mass_conservation_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [get_inversion_volume] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Consensus estimate optimisation with A factor: 10.0 and fs: 0
2023-03-07 12:35:08: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [mass_conservation_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [get_inversion_volume] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Consensus estimate optimisation with A factor: 9.927344372640782 and fs: 0
2023-03-07 12:35:08: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [mass_conservation_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [get_inversion_volume] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Consensus estimate optimisation with A factor: 9.503932371199765 and fs: 0
2023-03-07 12:35:08: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [mass_conservation_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [get_inversion_volume] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Consensus estimate optimisation with A factor: 9.551452033056764 and fs: 0
2023-03-07 12:35:08: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [mass_conservation_inversion] on 35 glaciers
2023-03-07 12:35:08: oggm.workflow: Execute entity tasks [get_inversion_volume] on 35 glaciers
2023-03-07 12:35:09: oggm.workflow: calibrate_inversion_from_consensus converged after 4 iterations and fs=0. The resulting Glen A factor is 9.503932371199765.
2023-03-07 12:35:09: oggm.workflow: Applying global task inversion_tasks on 35 glaciers
2023-03-07 12:35:09: oggm.workflow: Execute entity tasks [prepare_for_inversion] on 35 glaciers
2023-03-07 12:35:09: oggm.workflow: Execute entity tasks [mass_conservation_inversion] on 35 glaciers
2023-03-07 12:35:09: oggm.workflow: Execute entity tasks [get_inversion_volume] on 35 glaciers
cdf.sum()
vol_itmix_m3 6.423913e+07
vol_bsl_itmix_m3 0.000000e+00
vol_oggm_m3 6.425211e+07
dtype: float64
Note that here we calibrate the Glen A parameter to a value that is equal for all glaciers of gdirs (here \(A \sim 9.504\cdot A_0\)), i.e. we calibrate to match the total volume of all glaciers and not to match them individually.
cdf.iloc[:3]
| vol_itmix_m3 | vol_bsl_itmix_m3 | vol_oggm_m3 | |
|---|---|---|---|
| RGIId | |||
| RGI60-11.03208 | 1.564802e+07 | 0.0 | 1.972692e+07 |
| RGI60-11.03232 | 1.422478e+07 | 0.0 | 1.340392e+07 |
| RGI60-11.03222 | 5.158277e+06 | 0.0 | 6.604405e+06 |
just as a side note, “vol_bsl_itmix_m3” means volume below sea level and is therefore zero for these Alpine glaciers!
Distributed ice thickness#
The OGGM inversion and dynamical models use the “1D” flowline assumption: for some applications, you might want to use OGGM to create distributed ice thickness maps. Currently, OGGM implements two ways to “distribute” the flowline thicknesses, but only the simplest one works robustly:
# Distribute
workflow.execute_entity_task(tasks.distribute_thickness_per_altitude, gdirs);
2023-03-07 12:35:09: oggm.workflow: Execute entity tasks [distribute_thickness_per_altitude] on 35 glaciers
We just created a new output of the model, which we can access in the gridded_data file:
# xarray is an awesome library! Did you know about it?
import xarray as xr
import rioxarray as rioxr
ds = xr.open_dataset(gdirs[0].get_filepath('gridded_data'))
ds.distributed_thickness.plot();
Since some people find geotiff data easier to read than netCDF, OGGM also provides a tool to convert the variables in gridded_data.nc file to a geotiff file:
# save the distributed ice thickness into a geotiff file
workflow.execute_entity_task(tasks.gridded_data_var_to_geotiff, gdirs, varname='distributed_thickness')
# The default path of the geotiff file is in the glacier directory with the name "distributed_thickness.tif"
# Let's check if the file exists
for gdir in gdirs:
path = os.path.join(gdir.dir, 'distributed_thickness.tif')
assert os.path.exists(path)
2023-03-07 12:35:12: oggm.workflow: Execute entity tasks [gridded_data_var_to_geotiff] on 35 glaciers
# Open the last file with xarray's open_rasterio
rioxr.open_rasterio(path).plot();
In fact, tasks.gridded_data_var_to_geotiff() can save any variable in the gridded_data.nc file. The geotiff is named as the variable name with a .tif suffix. Have a try by yourself ;-)
Plot many glaciers on a map#
Let’s select a group of glaciers close to each other:
rgi_ids = ['RGI60-11.0{}'.format(i) for i in range(3205, 3211)]
sel_gdirs = [gdir for gdir in gdirs if gdir.rgi_id in rgi_ids]
graphics.plot_googlemap(sel_gdirs)
# you might need to install motionless if it is not yet in your environment
Using OGGM#
Since a recent PR (21.05.2020), OGGM can plot the thickness of a group of glaciers on a map:
graphics.plot_distributed_thickness(sel_gdirs)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
This is however not very useful because OGGM can only plot on a map as large as the local glacier map of the first glacier in the list. See this issue for a discussion about why.
Using salem#
Under the hood, OGGM uses salem to make the plots. Let’s do that for our case: it requires some manual tweaking, but it should be possible to automatize this better in the future.
Note: this also requires a very recent version of salem to work (21.05.2020)
import salem
# Make a grid covering the desired map extent
g = salem.mercator_grid(center_ll=(0.65, 42.64), extent=(4000, 4000))
# Create a map out of it
smap = salem.Map(g, countries=False)
# Add the glaciers outlines
for gdir in sel_gdirs:
crs = gdir.grid.center_grid
geom = gdir.read_pickle('geometries')
poly_pix = geom['polygon_pix']
smap.set_geometry(poly_pix, crs=crs, fc='none', zorder=2, linewidth=.2)
for l in poly_pix.interiors:
smap.set_geometry(l, crs=crs, color='black', linewidth=0.5)
f, ax = plt.subplots(figsize=(6, 6))
smap.visualize();
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
out = pickle.load(f)
# Now add the thickness data
for gdir in sel_gdirs:
grids_file = gdir.get_filepath('gridded_data')
with utils.ncDataset(grids_file) as nc:
vn = 'distributed_thickness'
thick = nc.variables[vn][:]
mask = nc.variables['glacier_mask'][:]
thick = np.where(mask, thick, np.NaN)
# The "overplot=True" is key here
# this needs a recent version of salem to run properly
smap.set_data(thick, crs=gdir.grid.center_grid, overplot=True)
# Set colorscale and other things
smap.set_cmap(graphics.OGGM_CMAPS['glacier_thickness'])
smap.set_plot_params(nlevels=256)
# Plot
f, ax = plt.subplots(figsize=(6, 6))
smap.visualize(ax=ax, cbar_title='Glacier thickness (m)');
