Ingest gridded products such as ice velocity into OGGM#

After running our OGGM experiments we often want to compare the model output to other gridded observations or maybe we want to use additional data sets that are not currently in the OGGM shop to calibrate parameters in the model (e.g. Glen A creep parameter, sliding parameter or the calving constant of proportionality). If you are looking on ways or ideas on how to do this, you are in the right tutorial!

In OGGM, a local map projection is defined for each glacier entity in the RGI inventory following the methods described in Maussion and others (2019). The model uses a Transverse Mercator projection centred on the glacier. A lot of data sets, especially those from Polar regions can have a different projections and if we are not careful, we would be making mistakes when we compare them with our model output or when we use such data sets to constrain our model experiments.

New in OGGM 1.6! We now offer preprocess directories where the data is available already. Visit OGGM as an accelerator for modelling and machine learning for more info.

First lets import the modules we need:

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import salem
from oggm import cfg, utils, workflow, tasks, graphics
cfg.initialize(logging_level='WARNING')
2024-02-02 17:05:15: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-02-02 17:05:15: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-02-02 17:05:15: oggm.cfg: Multiprocessing: using all available processors (N=4)
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-shop-on-Flowlines', reset=True)

Lets define the glaciers for the run#

rgi_ids = ['RGI60-14.06794']  # Baltoro
# The RGI version to use
# Size of the map around the glacier.
prepro_border = 80
# Degree of processing level. This is OGGM specific and for the shop 1 is the one you want
from_prepro_level = 3
# URL of the preprocessed gdirs
# we use elevation bands flowlines here
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5'
gdirs = workflow.init_glacier_directories(rgi_ids,
                                          from_prepro_level=from_prepro_level,
                                          prepro_base_url=base_url,
                                          prepro_border=prepro_border)
2024-02-02 17:05:15: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2024-02-02 17:05:15: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
gdir = gdirs[0]
graphics.plot_googlemap(gdir, figsize=(8, 7))
../../_images/dc39deb3bb6efdf057dae8b9f59a6147d00e9b78ba969036b4106f457121dffd.png

The gridded_data file in the glacier directory#

A lot of the data that the model use and produce for a glacier is stored under the glaciers directories in a NetCDF file called gridded_data:

fpath = gdir.get_filepath('gridded_data')
fpath
'/tmp/OGGM/OGGM-shop-on-Flowlines/per_glacier/RGI60-14/RGI60-14.06/RGI60-14.06794/gridded_data.nc'
with xr.open_dataset(fpath) as ds:
    ds = ds.load()
ds
<xarray.Dataset>
Dimensions:          (x: 479, y: 346)
Coordinates:
  * x                (x) float32 -4.764e+04 -4.744e+04 ... 4.776e+04 4.796e+04
  * y                (y) float32 3.992e+06 3.992e+06 ... 3.923e+06 3.923e+06
Data variables:
    topo             (y, x) float32 5.128e+03 4.999e+03 ... 5.291e+03 5.3e+03
    topo_smoothed    (y, x) float32 5.113e+03 5.019e+03 ... 5.292e+03 5.296e+03
    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
Attributes:
    author:         OGGM
    author_info:    Open Global Glacier Model
    pyproj_srs:     +proj=tmerc +lat_0=0 +lon_0=76.4047 +k=0.9996 +x_0=0 +y_0...
    max_h_dem:      8437.0
    min_h_dem:      2996.0
    max_h_glacier:  8437.0
    min_h_glacier:  3397.0
cmap=salem.get_cmap('topo')
ds.topo.plot(figsize=(7, 4), cmap=cmap);
../../_images/dc5b767a0ff93a0e29a4c32cb986e8b27e09ac9a4221ed433638c22e266dbb1b.png
ds.glacier_mask.plot();
../../_images/2050d334dda3f1e1279a13fac5c631eda1bfa9904c27b2819a19384fcdc24277.png

Add data from OGGM-Shop: bed topography data#

Additionally to the data produced by the model, the OGGM-Shop counts with routines that will automatically download and reproject other useful data sets into the glacier projection (For more information also check out this notebook). This data will be stored under the file described above.

OGGM can now 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-02-02 17:05:26: 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()

the cell below might take a while… be patient

# 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);
/usr/local/pyenv/versions/3.11.7/lib/python3.11/site-packages/salem/sio.py:429: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  dn = xarray_obj.dims.keys()
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/123d3d1b5d719b4a37e2f07cffd22c94ebfe39909538e4e2944a4cfbe361d600.png

OGGM-Shop: velocities#

We download data from Millan 2022 (see the shop).

If you want more velocity products, feel free to open a new topic on the OGGM issue tracker!

this will download severals large datasets depending on your connection, it might take some time

# attention downloads data!!!
from oggm.shop import millan22
workflow.execute_entity_task(millan22.velocity_to_gdir, gdirs);
2024-02-02 17:05:36: oggm.workflow: Execute entity tasks [velocity_to_gdir] on 1 glaciers

By applying the entity task velocity to gdir OGGM 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.

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

with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
    ds = ds.load()
ds
<xarray.Dataset>
Dimensions:                  (x: 479, y: 346)
Coordinates:
  * x                        (x) float32 -4.764e+04 -4.744e+04 ... 4.796e+04
  * y                        (y) float32 3.992e+06 3.992e+06 ... 3.923e+06
Data variables:
    topo                     (y, x) float32 5.128e+03 4.999e+03 ... 5.3e+03
    topo_smoothed            (y, x) float32 5.113e+03 5.019e+03 ... 5.296e+03
    topo_valid_mask          (y, x) int8 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
    glacier_ext              (y, x) int8 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
    consensus_ice_thickness  (y, x) float32 nan nan nan nan ... nan nan nan nan
    millan_v                 (y, x) float32 0.009646 0.0 8.727 ... 46.07 61.48
    millan_vx                (y, x) float32 nan nan 0.5949 ... -17.15 -28.61
    millan_vy                (y, x) float32 nan nan 8.707 ... 7.233 42.76 54.42
Attributes:
    author:         OGGM
    author_info:    Open Global Glacier Model
    pyproj_srs:     +proj=tmerc +lat_0=0 +lon_0=76.4047 +k=0.9996 +x_0=0 +y_0...
    max_h_dem:      8437.0
    min_h_dem:      2996.0
    max_h_glacier:  8437.0
    min_h_glacier:  3397.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);
/usr/local/pyenv/versions/3.11.7/lib/python3.11/site-packages/salem/sio.py:429: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  dn = xarray_obj.dims.keys()
ds
<xarray.Dataset>
Dimensions:                  (x: 479, y: 346)
Coordinates:
  * x                        (x) float32 -4.764e+04 -4.744e+04 ... 4.796e+04
  * y                        (y) float32 3.992e+06 3.992e+06 ... 3.923e+06
Data variables:
    topo                     (y, x) float32 5.128e+03 4.999e+03 ... 5.3e+03
    topo_smoothed            (y, x) float32 5.113e+03 5.019e+03 ... 5.296e+03
    topo_valid_mask          (y, x) int8 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
    glacier_ext              (y, x) int8 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
    consensus_ice_thickness  (y, x) float32 nan nan nan nan ... nan nan nan nan
    millan_v                 (y, x) float32 0.009646 0.0 8.727 ... 46.07 61.48
    millan_vx                (y, x) float32 nan nan 0.5949 ... -17.15 -28.61
    millan_vy                (y, x) float32 nan nan 8.707 ... 7.233 42.76 54.42
Attributes:
    author:         OGGM
    author_info:    Open Global Glacier Model
    pyproj_srs:     +proj=tmerc +lat_0=0 +lon_0=76.4047 +k=0.9996 +x_0=0 +y_0...
    max_h_dem:      8437.0
    min_h_dem:      2996.0
    max_h_glacier:  8437.0
    min_h_glacier:  3397.0
# get the velocity data
u = ds.millan_vx.where(ds.glacier_mask)
v = ds.millan_vy.where(ds.glacier_mask)
ws = (u**2 + v**2)**0.5

The .where(ds.glacier_mask) command will remove the data outside of the glacier outline.

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

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

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, 200, '200 m yr$^{-1}$',
                   labelpos='E', coordinates='axes')
../../_images/c225154e6a4c063610cf2c43c6a88491690ac3ca2bb0f6cee70e37eab4e95238.png

Bin the data in gridded_data into OGGM elevation bands#

Now that we have added new data to the gridded_data file, we can bin the data sets into the same elevation bands as OGGM by recomputing the elevation bands flowlines:

tasks.elevation_band_flowline(gdir,
                              bin_variables=['consensus_ice_thickness', 
                                             'millan_vx',
                                             'millan_vy'],
                              preserve_totals=[True, False, False]  # I"m actually not sure if preserving totals is meaningful with velocities - likely not
                              # NOTE: we could bin variables according to max() as well!
                              )

This created a csv in the glacier directory folder with the data binned to it:

import pandas as pd
df = pd.read_csv(gdir.get_filepath('elevation_band_flowline'), index_col=0)
df
area mean_elevation slope consensus_ice_thickness millan_vx millan_vy bin_elevation dx width
dis_along_flowline
60.171680 40000.0 8345.0800 0.244307 18.147762 4.305162 -1.510011 8355.0 120.343360 332.382274
140.402843 80000.0 8235.2370 0.642076 13.891147 2.600823 1.027054 8235.0 40.118965 1994.069381
173.939060 40000.0 8197.4900 0.838839 21.082249 10.003611 -6.313942 8205.0 26.953470 1484.038963
202.715851 40000.0 8150.4624 0.775496 25.654436 8.991845 3.440062 8145.0 30.600113 1307.184719
236.362152 40000.0 8075.4395 0.685384 14.702221 NaN NaN 8085.0 36.692490 1090.141345
... ... ... ... ... ... ... ... ... ...
42309.174751 760000.0 3525.2795 0.104130 139.481160 -4.504803 -6.100419 3525.0 287.058320 2647.545636
42584.145088 400000.0 3495.2273 0.113628 101.916240 -3.805607 -2.713950 3495.0 262.882350 1521.593187
42815.229424 200000.0 3472.7896 0.149415 98.673035 -0.934793 -5.369528 3465.0 199.286320 1003.581200
43010.617332 200000.0 3441.4329 0.155403 59.133278 -3.215284 -4.582561 3435.0 191.489500 1044.443679
43225.271721 120000.0 3407.1775 0.125483 32.996784 1.515356 -5.287875 3405.0 237.819270 504.584837

157 rows × 9 columns

df['obs_velocity'] = (df['millan_vx']**2 + df['millan_vy']**2)**0.5
df['bed_h'] = df['mean_elevation'] - df['consensus_ice_thickness']
df[['mean_elevation', 'bed_h', 'obs_velocity']].plot(figsize=(10, 4), secondary_y='obs_velocity');
../../_images/2740a2e34ca9c2a10f39f2cd4d06470191667c50e6ef6d3c43fc126c78351e18.png

The problem with this file is that it does not have a regular spacing. The numerical model needs regular spacing, which is why OGGM does this:

# This takes the csv file and prepares new 'inversion_flowlines.pkl' and created a new csv file with regular spacing
tasks.fixed_dx_elevation_band_flowline(gdir,
                                       bin_variables=['consensus_ice_thickness',
                                                      'millan_vx', 
                                                      'millan_vy'],
                                       preserve_totals=[True, False, False]
                                      )
df_regular = pd.read_csv(gdir.get_filepath('elevation_band_flowline', filesuffix='_fixed_dx'), index_col=0)
df_regular
widths_m area_m2 consensus_ice_thickness millan_vx millan_vy
dis_along_flowline
200.0 1350.505109 5.402020e+05 24.853092 9.087332 2.519514
600.0 1557.523130 6.230093e+05 12.058650 -2.843040 1.025657
1000.0 7529.727763 3.011891e+06 32.532274 4.897372 -1.923190
1400.0 15905.879873 6.362352e+06 29.213367 11.318087 -3.521698
1800.0 35987.372282 1.439495e+07 42.631259 0.019635 -0.396444
... ... ... ... ... ...
41400.0 2128.664211 8.514657e+05 188.035968 -7.261543 -24.703882
41800.0 1795.850103 7.183400e+05 163.141960 -9.733849 -13.058354
42200.0 2386.585691 9.546343e+05 148.052477 -5.063217 -7.519389
42600.0 1515.943567 6.063774e+05 100.202613 -3.608638 -2.896151
43000.0 1063.187354 4.252749e+05 60.383299 -3.091363 -4.625325

108 rows × 5 columns

The other variables have disappeared for saving space, but I think it would be nicer to have them here as well. We can grab them from the inversion flowlines (this could be better handled in OGGM):

fl = gdir.read_pickle('inversion_flowlines')[0]
df_regular['mean_elevation'] = fl.surface_h
df_regular['obs_velocity'] = (df_regular['millan_vx']**2 + df_regular['millan_vy']**2)**0.5
df_regular['consensus_bed_h'] = df_regular['mean_elevation'] - df_regular['consensus_ice_thickness']
df_regular[['mean_elevation', 'consensus_bed_h', 'obs_velocity']].plot(figsize=(10, 4),
                                                                       secondary_y='obs_velocity');
../../_images/f84d08712d4f3408700eb576b285d5ffe0fd3ef3c56b4d43517c13809aaca09d.png

OK so more or less the same but this time on a regular grid. Note that these now have the same length as the OGGM inversion flowlines, i.e. one can do stuff such as comparing what OGGM has done for the inversion:

inv = gdir.read_pickle('inversion_output')[0]
df_regular['OGGM_bed'] = df_regular['mean_elevation'] - inv['thick']
df_regular[['mean_elevation', 'consensus_bed_h', 'OGGM_bed']].plot();
../../_images/56cdff81dc2acce6967da047218baa7c1ac6767534ec9e55b393017df2813248.png

Compare velocities: at inversion#

OGGM already inverted the ice thickness based on some optimisation such as matching the regional totals. Which parameters did we use?

d = gdir.get_diagnostics()
d
{'dem_source': 'NASADEM',
 'flowline_type': 'elevation_band',
 'apparent_mb_from_any_mb_residual': 56.10000000000032,
 'inversion_glen_a': 6.42968271645714e-24,
 'inversion_fs': 0}
# OK set them so that we are consistent
cfg.PARAMS['inversion_glen_a'] = d['inversion_glen_a']
cfg.PARAMS['inversion_fs'] = d['inversion_fs']
2024-02-02 17:06:18: oggm.cfg: PARAMS['inversion_glen_a'] changed from `2.4e-24` to `6.42968271645714e-24`.
# Since we overwrote the inversion flowlines we have to do stuff again
tasks.mb_calibration_from_geodetic_mb(gdir,
                                      informed_threestep=True,
                                      overwrite_gdir=True)
tasks.apparent_mb_from_any_mb(gdir)
tasks.prepare_for_inversion(gdir)
tasks.mass_conservation_inversion(gdir)
tasks.compute_inversion_velocities(gdir)
inv = gdir.read_pickle('inversion_output')[0]
df_regular['OGGM_inversion_velocity'] = inv['u_surface']
df_regular[['obs_velocity', 'OGGM_inversion_velocity']].plot();
../../_images/b4738cc7c1772600167a51236182e92359801c620c3734c62501c36ed481ec84.png

The velocities are higher, because:

  • OGGM is at equilibrium

  • the comparison between the bulk velocities of OGGM and the gridded observed ones is not straightforward.

Velocities during the run#

TODO: here we could use velocities from the spinup run!!

Inversion velocities are for a glacier at equilibrium - this is not always meaningful. Lets do a run and store the velocities with time:

cfg.PARAMS['store_fl_diagnostics'] = True
2024-02-02 17:06:19: oggm.cfg: PARAMS['store_fl_diagnostics'] changed from `False` to `True`.
tasks.run_from_climate_data(gdir);
with xr.open_dataset(gdir.get_filepath('model_diagnostics')) as ds_diag:
    ds_diag = ds_diag.load()
ds_diag.volume_m3.plot();
../../_images/26e48d194d7d4ac3890c9dd733195d044131dca45fa827ee33cb1a5d24ad8422.png
with xr.open_dataset(gdir.get_filepath('fl_diagnostics'), group='fl_0') as ds_fl:
    ds_fl = ds_fl.load()
ds_fl
<xarray.Dataset>
Dimensions:              (dis_along_flowline: 179, time: 19)
Coordinates:
  * dis_along_flowline   (dis_along_flowline) float64 0.0 400.0 ... 7.12e+04
  * time                 (time) float64 2.002e+03 2.003e+03 ... 2.02e+03
Data variables: (12/13)
    point_lons           (dis_along_flowline) float64 75.88 75.88 ... 76.67
    point_lats           (dis_along_flowline) float64 36.07 36.07 ... 36.07
    bed_h                (dis_along_flowline) float64 8.131e+03 ... 3.007e+03
    volume_m3            (time, dis_along_flowline) float64 1.28e+07 ... 0.0
    volume_bsl_m3        (time, dis_along_flowline) float64 0.0 0.0 ... 0.0 0.0
    volume_bwl_m3        (time, dis_along_flowline) float64 0.0 0.0 ... 0.0 0.0
    ...                   ...
    thickness_m          (time, dis_along_flowline) float64 24.12 33.72 ... 0.0
    ice_velocity_myr     (time, dis_along_flowline) float64 nan nan ... 0.0 0.0
    calving_bucket_m3    (time) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    flux_divergence_myr  (time, dis_along_flowline) float64 nan nan ... 0.0 0.0
    climatic_mb_myr      (time, dis_along_flowline) float64 nan nan ... 0.0 0.0
    dhdt_myr             (time, dis_along_flowline) float64 nan nan ... 0.0 0.0
Attributes:
    class:                MixedBedFlowline
    map_dx:               200.0
    dx:                   2.0
    description:          OGGM model output
    oggm_version:         1.6.2.dev23+g5f8a8c5
    calendar:             365-day no leap
    creation_date:        2024-02-02 17:06:19
    water_level:          0
    glen_a:               6.42968271645714e-24
    fs:                   0
    mb_model_class:       MultipleFlowlineMassBalance
    mb_model_hemisphere:  nh
ds_fl.sel(time=[2003, 2010, 2020]).ice_velocity_myr.plot(hue='time');
../../_images/8501855edb675a7b9397345a70e00e4bc92b1996f1c19fa231424944613c7b3a.png
# The OGGM model flowlines also have the downstream lines
df_regular['OGGM_velocity_run_begin'] = ds_fl.sel(time=2003).ice_velocity_myr.data[:len(df_regular)]
df_regular['OGGM_velocity_run_end'] = ds_fl.sel(time=2020).ice_velocity_myr.data[:len(df_regular)]
df_regular[['obs_velocity', 'OGGM_inversion_velocity',
            'OGGM_velocity_run_begin', 'OGGM_velocity_run_end']].plot();
../../_images/e6e5705bc8927fe3fd39f5aa57bc6cd2107501a6c7169625573c2d4e05b26990.png

What’s next?#