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-08-25 21:14:37: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:14:37: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:14:37: 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-08-25 21:14:39: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2024-08-25 21:14:39: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers

gdir = gdirs[0]

graphics.plot_googlemap(gdir, figsize=(8, 7))

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> Size: 2MB
Dimensions:          (x: 479, y: 346)
Coordinates:
  * x                (x) float32 2kB -4.764e+04 -4.744e+04 ... 4.796e+04
  * y                (y) float32 1kB 3.992e+06 3.992e+06 ... 3.923e+06 3.923e+06
Data variables:
    topo             (y, x) float32 663kB 5.128e+03 4.999e+03 ... 5.3e+03
    topo_smoothed    (y, x) float32 663kB 5.113e+03 5.019e+03 ... 5.296e+03
    topo_valid_mask  (y, x) int8 166kB 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1
    glacier_mask     (y, x) int8 166kB 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0
    glacier_ext      (y, x) int8 166kB 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

xarray.Dataset

• Dimensions:
  • x: 479
  • y: 346
• Coordinates: (2)
  • x
    (x)
    float32
    -4.764e+04 -4.744e+04 ... 4.796e+04

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-47638.055, -47438.055, -47238.055, ...,  47561.945,  47761.945,
            47961.945], dtype=float32)

  • y
    (y)
    float32
    3.992e+06 3.992e+06 ... 3.923e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([3991927., 3991727., 3991527., ..., 3923327., 3923127., 3922927.],
          dtype=float32)

• Data variables: (5)
  • topo
    (y, x)
    float32
    5.128e+03 4.999e+03 ... 5.3e+03

    units :

    m

    long_name :

    DEM topography

    array([[5128., 4999., 4920., ..., 5500., 5563., 5641.],
           [5192., 5036., 4973., ..., 5423., 5449., 5520.],
           [5145., 5041., 5025., ..., 5260., 5313., 5396.],
           ...,
           [3985., 3874., 3787., ..., 5258., 5265., 5268.],
           [4032., 3927., 3822., ..., 5282., 5287., 5289.],
           [4113., 4021., 3934., ..., 5293., 5291., 5300.]], dtype=float32)

  • topo_smoothed
    (y, x)
    float32
    5.113e+03 5.019e+03 ... 5.296e+03

    units :

    m

    long_name :

    DEM topography smoothed with radius: 3e+02 m

    array([[5113.0107, 5019.37  , 4943.0776, ..., 5468.213 , 5543.3647,
            5599.1377],
           [5138.959 , 5046.7715, 4977.54  , ..., 5396.912 , 5451.8525,
            5503.462 ],
           [5128.378 , 5063.569 , 5022.3706, ..., 5261.8364, 5320.352 ,
            5378.3755],
           ...,
           [3953.9526, 3874.4246, 3787.1133, ..., 5262.384 , 5262.0576,
            5266.156 ],
           [4017.2664, 3936.9917, 3842.895 , ..., 5291.61  , 5283.133 ,
            5286.0977],
           [4075.7495, 4001.9124, 3910.3765, ..., 5304.3345, 5291.8545,
            5296.075 ]], dtype=float32)

  • topo_valid_mask
    (y, x)
    int8
    1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1

    units :

    -

    long_name :

    DEM validity mask according to geotiff input (1-0)

    array([[1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           ...,
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1]], dtype=int8)

  • glacier_mask
    (y, x)
    int8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    -

    long_name :

    Glacier mask

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           ...,
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)

  • glacier_ext
    (y, x)
    int8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    -

    long_name :

    Glacier external boundaries

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           ...,
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-47638.0546875, -47438.0546875, -47238.0546875, -47038.0546875,
           -46838.0546875, -46638.0546875, -46438.0546875, -46238.0546875,
           -46038.0546875, -45838.0546875,
           ...
            46161.9453125,  46361.9453125,  46561.9453125,  46761.9453125,
            46961.9453125,  47161.9453125,  47361.9453125,  47561.9453125,
            47761.9453125,  47961.9453125],
          dtype='float32', name='x', length=479))

  • y
    PandasIndex

    PandasIndex(Index([3991927.0, 3991727.0, 3991527.0, 3991327.0, 3991127.0, 3990927.0,
           3990727.0, 3990527.0, 3990327.0, 3990127.0,
           ...
           3924727.0, 3924527.0, 3924327.0, 3924127.0, 3923927.0, 3923727.0,
           3923527.0, 3923327.0, 3923127.0, 3922927.0],
          dtype='float32', name='y', length=346))

• Attributes: (7)

  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=0 +datum=WGS84 +units=m +no_defs

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

ds.glacier_mask.plot();

Merging the gridded_data files of multiple glacier directories

It is also possible to merge the gridded_data files of multiple glacier directories. Let’s try it with two glaciers:

rgi_ids_for_merge = ['RGI60-14.06794',  # Baltoro
                     'RGI60-14.07524',  # Siachen
                    ]
gdirs_for_merge = workflow.init_glacier_directories(rgi_ids_for_merge,
                                                    from_prepro_level=from_prepro_level,
                                                    prepro_base_url=base_url,
                                                    prepro_border=prepro_border)

2024-08-25 21:14:50: oggm.workflow: init_glacier_directories from prepro level 3 on 2 glaciers.
2024-08-25 21:14:50: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 2 glaciers

ds_merged = workflow.merge_gridded_data(
    gdirs_for_merge,
    output_folder=None,  # by default the final file is saved at cfg.PATHS['working_dir']
    output_filename='gridded_data_merged',  # the default file is saved as gridded_data_merged.nc
    included_variables='all',  # you also can provide a list of variables here
    add_topography=False,  # here we can add topography for the new extend
    reset=False,  # set to True if you want to overwrite an already existing file (for playing around)
)

2024-08-25 21:14:50: oggm.workflow: Applying global task merge_gridded_data on 2 glaciers
2024-08-25 21:14:50: oggm.workflow: Execute entity tasks [reproject_gridded_data_variable_to_grid] on 2 glaciers
2024-08-25 21:14:50: oggm.workflow: Execute entity tasks [reproject_gridded_data_variable_to_grid] on 2 glaciers

ds_merged

<xarray.Dataset> Size: 3MB
Dimensions:       (x: 744, y: 582)
Coordinates:
  * x             (x) float32 3kB -4.774e+04 -4.754e+04 ... 1.007e+05 1.009e+05
  * y             (y) float32 2kB 3.992e+06 3.992e+06 ... 3.876e+06 3.876e+06
Data variables:
    glacier_mask  (y, x) float32 2MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    glacier_ext   (y, x) float32 2MB 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_...
    nr_of_merged_glaciers:  2
    rgi_ids:                ['RGI60-14.06794', 'RGI60-14.07524']

xarray.Dataset

• Dimensions:
  • x: 744
  • y: 582
• Coordinates: (2)
  • x
    (x)
    float32
    -4.774e+04 -4.754e+04 ... 1.009e+05

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-47738.055, -47538.055, -47338.055, ..., 100461.945, 100661.945,
           100861.945], dtype=float32)

  • y
    (y)
    float32
    3.992e+06 3.992e+06 ... 3.876e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([3992027., 3991827., 3991627., ..., 3876227., 3876027., 3875827.],
          dtype=float32)

• Data variables: (2)
  • glacier_mask
    (y, x)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    -

    long_name :

    Glacier mask

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)

  • glacier_ext
    (y, x)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    units :

    -

    long_name :

    Glacier external boundaries

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-47738.0546875, -47538.0546875, -47338.0546875, -47138.0546875,
           -46938.0546875, -46738.0546875, -46538.0546875, -46338.0546875,
           -46138.0546875, -45938.0546875,
           ...
            99061.9453125,  99261.9453125,  99461.9453125,  99661.9453125,
            99861.9453125, 100061.9453125, 100261.9453125, 100461.9453125,
           100661.9453125, 100861.9453125],
          dtype='float32', name='x', length=744))

  • y
    PandasIndex

    PandasIndex(Index([3992027.0, 3991827.0, 3991627.0, 3991427.0, 3991227.0, 3991027.0,
           3990827.0, 3990627.0, 3990427.0, 3990227.0,
           ...
           3877627.0, 3877427.0, 3877227.0, 3877027.0, 3876827.0, 3876627.0,
           3876427.0, 3876227.0, 3876027.0, 3875827.0],
          dtype='float32', name='y', length=582))

• Attributes: (5)

  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=0 +datum=WGS84 +units=m +no_defs

  nr_of_merged_glaciers :

  2

  rgi_ids :

  ['RGI60-14.06794', 'RGI60-14.07524']

ds_merged.glacier_mask.plot()

<matplotlib.collections.QuadMesh at 0x7fee99d339d0>

As you can see, topography is not included because this process requires some time. The reason is that topography is not simply reprojected from the existing gridded_data; instead, it is redownloaded from the source and processed for the new ‘merged extent’. You can include topography by setting add_topography=True or specifying the DEM source add_topography='dem_source'.

There are additional options available, such as use_glacier_mask or preserve_totals (for preserving the total volume when merging distributed thickness data). For a more comprehensive explanation of these options, refer to the Docstring (by uncommenting the line below).

# for more available options uncomment the following line
# workflow.merge_gridded_data?

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-08-25 21:14:51: 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);

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: 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-08-25 21:15:02: 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> Size: 4MB
Dimensions:                  (x: 479, y: 346)
Coordinates:
  * x                        (x) float32 2kB -4.764e+04 -4.744e+04 ... 4.796e+04
  * y                        (y) float32 1kB 3.992e+06 3.992e+06 ... 3.923e+06
Data variables:
    topo                     (y, x) float32 663kB 5.128e+03 ... 5.3e+03
    topo_smoothed            (y, x) float32 663kB 5.113e+03 ... 5.296e+03
    topo_valid_mask          (y, x) int8 166kB 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1
    glacier_mask             (y, x) int8 166kB 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
    glacier_ext              (y, x) int8 166kB 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
    consensus_ice_thickness  (y, x) float32 663kB nan nan nan ... nan nan nan
    millan_v                 (y, x) float32 663kB 0.009646 0.0 ... 46.07 61.48
    millan_vx                (y, x) float32 663kB nan nan ... -17.15 -28.61
    millan_vy                (y, x) float32 663kB nan nan 8.707 ... 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

xarray.Dataset

• Dimensions:
  • x: 479
  • y: 346
• Coordinates: (2)
  • x
    (x)
    float32
    -4.764e+04 -4.744e+04 ... 4.796e+04

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-47638.055, -47438.055, -47238.055, ...,  47561.945,  47761.945,
            47961.945], dtype=float32)

  • y
    (y)
    float32
    3.992e+06 3.992e+06 ... 3.923e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([3991927., 3991727., 3991527., ..., 3923327., 3923127., 3922927.],
          dtype=float32)

• Data variables: (9)
  • topo
    (y, x)
    float32
    5.128e+03 4.999e+03 ... 5.3e+03

    units :

    m

    long_name :

    DEM topography

    array([[5128., 4999., 4920., ..., 5500., 5563., 5641.],
           [5192., 5036., 4973., ..., 5423., 5449., 5520.],
           [5145., 5041., 5025., ..., 5260., 5313., 5396.],
           ...,
           [3985., 3874., 3787., ..., 5258., 5265., 5268.],
           [4032., 3927., 3822., ..., 5282., 5287., 5289.],
           [4113., 4021., 3934., ..., 5293., 5291., 5300.]], dtype=float32)

  • topo_smoothed
    (y, x)
    float32
    5.113e+03 5.019e+03 ... 5.296e+03

    units :

    m

    long_name :

    DEM topography smoothed with radius: 3e+02 m

    array([[5113.0107, 5019.37  , 4943.0776, ..., 5468.213 , 5543.3647,
            5599.1377],
           [5138.959 , 5046.7715, 4977.54  , ..., 5396.912 , 5451.8525,
            5503.462 ],
           [5128.378 , 5063.569 , 5022.3706, ..., 5261.8364, 5320.352 ,
            5378.3755],
           ...,
           [3953.9526, 3874.4246, 3787.1133, ..., 5262.384 , 5262.0576,
            5266.156 ],
           [4017.2664, 3936.9917, 3842.895 , ..., 5291.61  , 5283.133 ,
            5286.0977],
           [4075.7495, 4001.9124, 3910.3765, ..., 5304.3345, 5291.8545,
            5296.075 ]], dtype=float32)

  • topo_valid_mask
    (y, x)
    int8
    1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1

    units :

    -

    long_name :

    DEM validity mask according to geotiff input (1-0)

    array([[1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           ...,
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1]], dtype=int8)

  • glacier_mask
    (y, x)
    int8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    -

    long_name :

    Glacier mask

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           ...,
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)

  • glacier_ext
    (y, x)
    int8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    -

    long_name :

    Glacier external boundaries

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           ...,
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)

  • consensus_ice_thickness
    (y, x)
    float32
    nan nan nan nan ... nan nan nan nan

    units :

    m

    long_name :

    Ice thickness from the consensus estimate

    base_url :

    https://cluster.klima.uni-bremen.de/~fmaussion/icevol/composite/

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

  • millan_v
    (y, x)
    float32
    0.009646 0.0 8.727 ... 46.07 61.48

    long_name :

    Ice velocity from Millan et al. 2022

    data_source :

    RGI-13-15.6_2022February09

    units :

    m

    array([[9.6464921e-03, 0.0000000e+00, 8.7274847e+00, ..., 0.0000000e+00,
            0.0000000e+00, 0.0000000e+00],
           [0.0000000e+00, 3.3396182e+00, 1.9843628e+01, ..., 0.0000000e+00,
            0.0000000e+00, 0.0000000e+00],
           [0.0000000e+00, 4.4036991e+01, 1.9893900e+01, ..., 0.0000000e+00,
            0.0000000e+00, 0.0000000e+00],
           ...,
           [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ..., 4.4023067e+01,
            5.2307713e+01, 5.8931698e+01],
           [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ..., 3.1996721e+01,
            5.2488129e+01, 5.7668842e+01],
           [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ..., 1.2152705e+01,
            4.6074623e+01, 6.1482475e+01]], dtype=float32)

  • millan_vx
    (y, x)
    float32
    nan nan 0.5949 ... -17.15 -28.61

    units :

    m

    long_name :

    Ice velocity in map x direction from Millan et al. 2022

    data_source :

    RGI-13-15.6_2022February09

    array([[        nan,         nan,   0.5949099, ...,         nan,
                    nan,         nan],
           [        nan,         nan,  11.420494 , ...,         nan,
                    nan,         nan],
           [        nan,  43.813683 ,   4.029796 , ...,         nan,
                    nan,         nan],
           ...,
           [        nan,         nan,         nan, ..., -11.01701  ,
             -5.947851 , -21.181208 ],
           [        nan,         nan,         nan, ...,  -9.892379 ,
            -17.386967 , -22.882816 ],
           [        nan,         nan,         nan, ...,  -9.765955 ,
            -17.154419 , -28.61151  ]], dtype=float32)

  • millan_vy
    (y, x)
    float32
    nan nan 8.707 ... 7.233 42.76 54.42

    units :

    m

    long_name :

    Ice velocity in map y direction from Millan et al. 2022

    data_source :

    RGI-13-15.6_2022February09

    array([[       nan,        nan,  8.707185 , ...,        nan,        nan,
                   nan],
           [       nan,        nan, 16.227812 , ...,        nan,        nan,
                   nan],
           [       nan, -4.4291787, 19.481478 , ...,        nan,        nan,
                   nan],
           ...,
           [       nan,        nan,        nan, ..., 42.622246 , 51.968452 ,
            54.99365  ],
           [       nan,        nan,        nan, ..., 30.429115 , 49.52471  ,
            52.934605 ],
           [       nan,        nan,        nan, ...,  7.2328672, 42.762093 ,
            54.41945  ]], dtype=float32)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-47638.0546875, -47438.0546875, -47238.0546875, -47038.0546875,
           -46838.0546875, -46638.0546875, -46438.0546875, -46238.0546875,
           -46038.0546875, -45838.0546875,
           ...
            46161.9453125,  46361.9453125,  46561.9453125,  46761.9453125,
            46961.9453125,  47161.9453125,  47361.9453125,  47561.9453125,
            47761.9453125,  47961.9453125],
          dtype='float32', name='x', length=479))

  • y
    PandasIndex

    PandasIndex(Index([3991927.0, 3991727.0, 3991527.0, 3991327.0, 3991127.0, 3990927.0,
           3990727.0, 3990527.0, 3990327.0, 3990127.0,
           ...
           3924727.0, 3924527.0, 3924327.0, 3924127.0, 3923927.0, 3923727.0,
           3923527.0, 3923327.0, 3923127.0, 3922927.0],
          dtype='float32', name='y', length=346))

• Attributes: (7)

  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=0 +datum=WGS84 +units=m +no_defs

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

ds

<xarray.Dataset> Size: 4MB
Dimensions:                  (x: 479, y: 346)
Coordinates:
  * x                        (x) float32 2kB -4.764e+04 -4.744e+04 ... 4.796e+04
  * y                        (y) float32 1kB 3.992e+06 3.992e+06 ... 3.923e+06
Data variables:
    topo                     (y, x) float32 663kB 5.128e+03 ... 5.3e+03
    topo_smoothed            (y, x) float32 663kB 5.113e+03 ... 5.296e+03
    topo_valid_mask          (y, x) int8 166kB 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1
    glacier_mask             (y, x) int8 166kB 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
    glacier_ext              (y, x) int8 166kB 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0
    consensus_ice_thickness  (y, x) float32 663kB nan nan nan ... nan nan nan
    millan_v                 (y, x) float32 663kB 0.009646 0.0 ... 46.07 61.48
    millan_vx                (y, x) float32 663kB nan nan ... -17.15 -28.61
    millan_vy                (y, x) float32 663kB nan nan 8.707 ... 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

xarray.Dataset

• Dimensions:
  • x: 479
  • y: 346
• Coordinates: (2)
  • x
    (x)
    float32
    -4.764e+04 -4.744e+04 ... 4.796e+04

    units :

    m

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    array([-47638.055, -47438.055, -47238.055, ...,  47561.945,  47761.945,
            47961.945], dtype=float32)

  • y
    (y)
    float32
    3.992e+06 3.992e+06 ... 3.923e+06

    units :

    m

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    array([3991927., 3991727., 3991527., ..., 3923327., 3923127., 3922927.],
          dtype=float32)

• Data variables: (9)
  • topo
    (y, x)
    float32
    5.128e+03 4.999e+03 ... 5.3e+03

    units :

    m

    long_name :

    DEM topography

    array([[5128., 4999., 4920., ..., 5500., 5563., 5641.],
           [5192., 5036., 4973., ..., 5423., 5449., 5520.],
           [5145., 5041., 5025., ..., 5260., 5313., 5396.],
           ...,
           [3985., 3874., 3787., ..., 5258., 5265., 5268.],
           [4032., 3927., 3822., ..., 5282., 5287., 5289.],
           [4113., 4021., 3934., ..., 5293., 5291., 5300.]], dtype=float32)

  • topo_smoothed
    (y, x)
    float32
    5.113e+03 5.019e+03 ... 5.296e+03

    units :

    m

    long_name :

    DEM topography smoothed with radius: 3e+02 m

    array([[5113.0107, 5019.37  , 4943.0776, ..., 5468.213 , 5543.3647,
            5599.1377],
           [5138.959 , 5046.7715, 4977.54  , ..., 5396.912 , 5451.8525,
            5503.462 ],
           [5128.378 , 5063.569 , 5022.3706, ..., 5261.8364, 5320.352 ,
            5378.3755],
           ...,
           [3953.9526, 3874.4246, 3787.1133, ..., 5262.384 , 5262.0576,
            5266.156 ],
           [4017.2664, 3936.9917, 3842.895 , ..., 5291.61  , 5283.133 ,
            5286.0977],
           [4075.7495, 4001.9124, 3910.3765, ..., 5304.3345, 5291.8545,
            5296.075 ]], dtype=float32)

  • topo_valid_mask
    (y, x)
    int8
    1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1

    units :

    -

    long_name :

    DEM validity mask according to geotiff input (1-0)

    array([[1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           ...,
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1],
           [1, 1, 1, ..., 1, 1, 1]], dtype=int8)

  • glacier_mask
    (y, x)
    int8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    -

    long_name :

    Glacier mask

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           ...,
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)

  • glacier_ext
    (y, x)
    int8
    0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0

    units :

    -

    long_name :

    Glacier external boundaries

    array([[0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           ...,
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0],
           [0, 0, 0, ..., 0, 0, 0]], dtype=int8)

  • consensus_ice_thickness
    (y, x)
    float32
    nan nan nan nan ... nan nan nan nan

    units :

    m

    long_name :

    Ice thickness from the consensus estimate

    base_url :

    https://cluster.klima.uni-bremen.de/~fmaussion/icevol/composite/

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

  • millan_v
    (y, x)
    float32
    0.009646 0.0 8.727 ... 46.07 61.48

    long_name :

    Ice velocity from Millan et al. 2022

    data_source :

    RGI-13-15.6_2022February09

    units :

    m

    array([[9.6464921e-03, 0.0000000e+00, 8.7274847e+00, ..., 0.0000000e+00,
            0.0000000e+00, 0.0000000e+00],
           [0.0000000e+00, 3.3396182e+00, 1.9843628e+01, ..., 0.0000000e+00,
            0.0000000e+00, 0.0000000e+00],
           [0.0000000e+00, 4.4036991e+01, 1.9893900e+01, ..., 0.0000000e+00,
            0.0000000e+00, 0.0000000e+00],
           ...,
           [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ..., 4.4023067e+01,
            5.2307713e+01, 5.8931698e+01],
           [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ..., 3.1996721e+01,
            5.2488129e+01, 5.7668842e+01],
           [0.0000000e+00, 0.0000000e+00, 0.0000000e+00, ..., 1.2152705e+01,
            4.6074623e+01, 6.1482475e+01]], dtype=float32)

  • millan_vx
    (y, x)
    float32
    nan nan 0.5949 ... -17.15 -28.61

    units :

    m

    long_name :

    Ice velocity in map x direction from Millan et al. 2022

    data_source :

    RGI-13-15.6_2022February09

    array([[        nan,         nan,   0.5949099, ...,         nan,
                    nan,         nan],
           [        nan,         nan,  11.420494 , ...,         nan,
                    nan,         nan],
           [        nan,  43.813683 ,   4.029796 , ...,         nan,
                    nan,         nan],
           ...,
           [        nan,         nan,         nan, ..., -11.01701  ,
             -5.947851 , -21.181208 ],
           [        nan,         nan,         nan, ...,  -9.892379 ,
            -17.386967 , -22.882816 ],
           [        nan,         nan,         nan, ...,  -9.765955 ,
            -17.154419 , -28.61151  ]], dtype=float32)

  • millan_vy
    (y, x)
    float32
    nan nan 8.707 ... 7.233 42.76 54.42

    units :

    m

    long_name :

    Ice velocity in map y direction from Millan et al. 2022

    data_source :

    RGI-13-15.6_2022February09

    array([[       nan,        nan,  8.707185 , ...,        nan,        nan,
                   nan],
           [       nan,        nan, 16.227812 , ...,        nan,        nan,
                   nan],
           [       nan, -4.4291787, 19.481478 , ...,        nan,        nan,
                   nan],
           ...,
           [       nan,        nan,        nan, ..., 42.622246 , 51.968452 ,
            54.99365  ],
           [       nan,        nan,        nan, ..., 30.429115 , 49.52471  ,
            52.934605 ],
           [       nan,        nan,        nan, ...,  7.2328672, 42.762093 ,
            54.41945  ]], dtype=float32)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-47638.0546875, -47438.0546875, -47238.0546875, -47038.0546875,
           -46838.0546875, -46638.0546875, -46438.0546875, -46238.0546875,
           -46038.0546875, -45838.0546875,
           ...
            46161.9453125,  46361.9453125,  46561.9453125,  46761.9453125,
            46961.9453125,  47161.9453125,  47361.9453125,  47561.9453125,
            47761.9453125,  47961.9453125],
          dtype='float32', name='x', length=479))

  • y
    PandasIndex

    PandasIndex(Index([3991927.0, 3991727.0, 3991527.0, 3991327.0, 3991127.0, 3990927.0,
           3990727.0, 3990527.0, 3990327.0, 3990127.0,
           ...
           3924727.0, 3924527.0, 3924327.0, 3924127.0, 3923927.0, 3923727.0,
           3923527.0, 3923327.0, 3923127.0, 3922927.0],
          dtype='float32', name='y', length=346))

• Attributes: (7)

  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=0 +datum=WGS84 +units=m +no_defs

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

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

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
0.0	1350.505109	5.402020e+05	24.853092	9.087332	2.519514
400.0	1557.523130	6.230093e+05	12.058650	-2.843040	1.025657
800.0	7529.727763	3.011891e+06	32.532274	4.897372	-1.923190
1200.0	15905.879873	6.362352e+06	29.213367	11.318087	-3.521698
1600.0	35987.372282	1.439495e+07	42.631259	0.019635	-0.396444
...	...	...	...	...	...
41200.0	2128.664211	8.514657e+05	188.035968	-7.261543	-24.703882
41600.0	1795.850103	7.183400e+05	163.141960	-9.733849	-13.058354
42000.0	2386.585691	9.546343e+05	148.052477	-5.063217	-7.519389
42400.0	1515.943567	6.063774e+05	100.202613	-3.608638	-2.896151
42800.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');

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

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-08-25 21:16:00: 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();

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-08-25 21:16:01: 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();

with xr.open_dataset(gdir.get_filepath('fl_diagnostics'), group='fl_0') as ds_fl:
    ds_fl = ds_fl.load()

ds_fl

<xarray.Dataset> Size: 251kB
Dimensions:              (dis_along_flowline: 179, time: 19)
Coordinates:
  * dis_along_flowline   (dis_along_flowline) float64 1kB 0.0 400.0 ... 7.12e+04
  * time                 (time) float64 152B 2.002e+03 2.003e+03 ... 2.02e+03
Data variables: (12/13)
    point_lons           (dis_along_flowline) float64 1kB 75.88 75.88 ... 76.67
    point_lats           (dis_along_flowline) float64 1kB 36.07 36.07 ... 36.07
    bed_h                (dis_along_flowline) float64 1kB 8.131e+03 ... 3.007...
    volume_m3            (time, dis_along_flowline) float64 27kB 1.28e+07 ......
    volume_bsl_m3        (time, dis_along_flowline) float64 27kB 0.0 0.0 ... 0.0
    volume_bwl_m3        (time, dis_along_flowline) float64 27kB 0.0 0.0 ... 0.0
    ...                   ...
    thickness_m          (time, dis_along_flowline) float64 27kB 24.12 ... 0.0
    ice_velocity_myr     (time, dis_along_flowline) float64 27kB nan nan ... 0.0
    calving_bucket_m3    (time) float64 152B 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    flux_divergence_myr  (time, dis_along_flowline) float64 27kB nan nan ... 0.0
    climatic_mb_myr      (time, dis_along_flowline) float64 27kB nan nan ... 0.0
    dhdt_myr             (time, dis_along_flowline) float64 27kB nan nan ... 0.0
Attributes:
    calendar:             365-day no leap
    water_level:          0
    dx:                   2.0
    creation_date:        2024-08-25 21:16:01
    fs:                   0
    mb_model_class:       MultipleFlowlineMassBalance
    map_dx:               200.0
    class:                MixedBedFlowline
    glen_a:               6.42968271645714e-24
    oggm_version:         1.6.3.dev3+g865475d
    mb_model_hemisphere:  nh
    description:          OGGM model output

xarray.Dataset

• Dimensions:
  • dis_along_flowline: 179
  • time: 19
• Coordinates: (2)
  • dis_along_flowline
    (dis_along_flowline)
    float64
    0.0 400.0 ... 7.08e+04 7.12e+04

    array([    0.,   400.,   800.,  1200.,  1600.,  2000.,  2400.,  2800.,  3200.,
            3600.,  4000.,  4400.,  4800.,  5200.,  5600.,  6000.,  6400.,  6800.,
            7200.,  7600.,  8000.,  8400.,  8800.,  9200.,  9600., 10000., 10400.,
           10800., 11200., 11600., 12000., 12400., 12800., 13200., 13600., 14000.,
           14400., 14800., 15200., 15600., 16000., 16400., 16800., 17200., 17600.,
           18000., 18400., 18800., 19200., 19600., 20000., 20400., 20800., 21200.,
           21600., 22000., 22400., 22800., 23200., 23600., 24000., 24400., 24800.,
           25200., 25600., 26000., 26400., 26800., 27200., 27600., 28000., 28400.,
           28800., 29200., 29600., 30000., 30400., 30800., 31200., 31600., 32000.,
           32400., 32800., 33200., 33600., 34000., 34400., 34800., 35200., 35600.,
           36000., 36400., 36800., 37200., 37600., 38000., 38400., 38800., 39200.,
           39600., 40000., 40400., 40800., 41200., 41600., 42000., 42400., 42800.,
           43200., 43600., 44000., 44400., 44800., 45200., 45600., 46000., 46400.,
           46800., 47200., 47600., 48000., 48400., 48800., 49200., 49600., 50000.,
           50400., 50800., 51200., 51600., 52000., 52400., 52800., 53200., 53600.,
           54000., 54400., 54800., 55200., 55600., 56000., 56400., 56800., 57200.,
           57600., 58000., 58400., 58800., 59200., 59600., 60000., 60400., 60800.,
           61200., 61600., 62000., 62400., 62800., 63200., 63600., 64000., 64400.,
           64800., 65200., 65600., 66000., 66400., 66800., 67200., 67600., 68000.,
           68400., 68800., 69200., 69600., 70000., 70400., 70800., 71200.])

  • time
    (time)
    float64
    2.002e+03 2.003e+03 ... 2.02e+03

    description :

    Floating hydrological year

    array([2002., 2003., 2004., 2005., 2006., 2007., 2008., 2009., 2010., 2011.,
           2012., 2013., 2014., 2015., 2016., 2017., 2018., 2019., 2020.])

• Data variables: (13)
  • point_lons
    (dis_along_flowline)
    float64
    75.88 75.88 75.88 ... 76.66 76.67

    description :

    Longitude along the flowline

    unit :

    deg

    array([75.87566708, 75.88010902, 75.88455096, 75.8889929 , 75.89343485,
           75.8978768 , 75.90231876, 75.90676072, 75.91120268, 75.91564465,
           75.92008662, 75.9245286 , 75.92897058, 75.93341256, 75.93785455,
           75.94229654, 75.94673854, 75.95118054, 75.95562254, 75.96006455,
           75.96450656, 75.96894857, 75.97339059, 75.97783261, 75.98227463,
           75.98671666, 75.99115869, 75.99560073, 76.00004276, 76.0044848 ,
           76.00892685, 76.0133689 , 76.01781095, 76.022253  , 76.02669506,
           76.03113712, 76.03557918, 76.04002125, 76.04446331, 76.04890539,
           76.05334746, 76.05778954, 76.06223162, 76.0666737 , 76.07111579,
           76.07555788, 76.07999997, 76.08444206, 76.08888416, 76.09332626,
           76.09776836, 76.10221046, 76.10665257, 76.11109468, 76.11553679,
           76.1199789 , 76.12442102, 76.12886314, 76.13330526, 76.13774738,
           76.1421895 , 76.14663163, 76.15107376, 76.15551589, 76.15995802,
           76.16440016, 76.1688423 , 76.17328444, 76.17772658, 76.18216872,
           76.18661086, 76.19105301, 76.19549516, 76.19993731, 76.20437946,
           76.20882161, 76.21326377, 76.21770592, 76.22214808, 76.22659024,
           76.2310324 , 76.23547456, 76.23991672, 76.24435889, 76.24880106,
           76.25324322, 76.25768539, 76.26212756, 76.26656973, 76.2710119 ,
           76.27545408, 76.27989625, 76.28433843, 76.2887806 , 76.29322278,
           76.29766496, 76.30210713, 76.30654931, 76.31099149, 76.31543368,
           76.31987586, 76.32431804, 76.32876022, 76.33320241, 76.33764459,
           76.34208678, 76.34652896, 76.35097115, 76.35541333, 76.35985552,
           76.36429771, 76.3687399 , 76.37318208, 76.37762427, 76.38206646,
           76.38650865, 76.39095084, 76.39539303, 76.39983521, 76.4042774 ,
           76.40871959, 76.41316178, 76.41760397, 76.42204616, 76.42648835,
           76.43093053, 76.43537272, 76.43981491, 76.4442571 , 76.44869928,
           76.45314147, 76.45758366, 76.46202584, 76.46646803, 76.47091021,
           76.4753524 , 76.47979458, 76.48423677, 76.48867895, 76.49312113,
           76.49756331, 76.50200549, 76.50644767, 76.51088985, 76.51533203,
           76.51977421, 76.52421638, 76.52865856, 76.53310073, 76.53754291,
           76.54198508, 76.54642725, 76.55086942, 76.55531159, 76.55975375,
           76.56419592, 76.56863809, 76.57308025, 76.57752241, 76.58196457,
           76.58640673, 76.59084889, 76.59529105, 76.5997332 , 76.60417535,
           76.60861751, 76.61305966, 76.6175018 , 76.62194395, 76.6263861 ,
           76.63082824, 76.63527038, 76.63971252, 76.64415466, 76.64859679,
           76.65303893, 76.65748106, 76.66192319, 76.66636531])

  • point_lats
    (dis_along_flowline)
    float64
    36.07 36.07 36.07 ... 36.07 36.07

    description :

    Latitude along the flowline

    unit :

    deg

    array([36.07076715, 36.07078668, 36.07080604, 36.07082523, 36.07084426,
           36.07086312, 36.07088182, 36.07090036, 36.07091873, 36.07093694,
           36.07095498, 36.07097285, 36.07099057, 36.07100812, 36.0710255 ,
           36.07104272, 36.07105977, 36.07107666, 36.07109339, 36.07110995,
           36.07112634, 36.07114257, 36.07115864, 36.07117454, 36.07119028,
           36.07120585, 36.07122126, 36.0712365 , 36.07125158, 36.07126649,
           36.07128124, 36.07129583, 36.07131025, 36.07132451, 36.0713386 ,
           36.07135252, 36.07136628, 36.07137988, 36.07139331, 36.07140658,
           36.07141969, 36.07143262, 36.0714454 , 36.07145801, 36.07147045,
           36.07148273, 36.07149485, 36.0715068 , 36.07151859, 36.07153021,
           36.07154167, 36.07155296, 36.07156409, 36.07157505, 36.07158585,
           36.07159648, 36.07160695, 36.07161726, 36.0716274 , 36.07163737,
           36.07164718, 36.07165683, 36.07166631, 36.07167563, 36.07168478,
           36.07169377, 36.07170259, 36.07171125, 36.07171974, 36.07172807,
           36.07173624, 36.07174424, 36.07175207, 36.07175975, 36.07176725,
           36.07177459, 36.07178177, 36.07178878, 36.07179563, 36.07180231,
           36.07180883, 36.07181519, 36.07182138, 36.0718274 , 36.07183326,
           36.07183896, 36.07184449, 36.07184985, 36.07185505, 36.07186009,
           36.07186496, 36.07186967, 36.07187421, 36.07187859, 36.07188281,
           36.07188685, 36.07189074, 36.07189446, 36.07189801, 36.0719014 ,
           36.07190463, 36.07190769, 36.07191059, 36.07191332, 36.07191589,
           36.07191829, 36.07192053, 36.0719226 , 36.07192451, 36.07192626,
           36.07192783, 36.07192925, 36.0719305 , 36.07193159, 36.07193251,
           36.07193326, 36.07193386, 36.07193428, 36.07193455, 36.07193464,
           36.07193458, 36.07193435, 36.07193395, 36.07193339, 36.07193266,
           36.07193177, 36.07193072, 36.0719295 , 36.07192812, 36.07192657,
           36.07192486, 36.07192298, 36.07192094, 36.07191873, 36.07191636,
           36.07191382, 36.07191112, 36.07190826, 36.07190523, 36.07190203,
           36.07189867, 36.07189515, 36.07189146, 36.07188761, 36.07188359,
           36.07187941, 36.07187506, 36.07187055, 36.07186587, 36.07186103,
           36.07185602, 36.07185085, 36.07184552, 36.07184002, 36.07183436,
           36.07182853, 36.07182253, 36.07181638, 36.07181005, 36.07180357,
           36.07179691, 36.0717901 , 36.07178312, 36.07177597, 36.07176866,
           36.07176119, 36.07175355, 36.07174574, 36.07173777, 36.07172964,
           36.07172134, 36.07171288, 36.07170425, 36.07169546, 36.0716865 ,
           36.07167738, 36.0716681 , 36.07165865, 36.07164903])

  • bed_h
    (dis_along_flowline)
    float64
    8.131e+03 7.694e+03 ... 3.007e+03

    description :

    Bed elevation along the flowline

    unit :

    m

    array([8130.776148  , 7693.73914873, 7513.53649227, 7249.84999508,
           7045.74897083, 6799.45212016, 6561.62473351, 6361.77759514,
           6101.12475002, 5881.21284792, 5685.73069301, 5520.6492943 ,
           5401.15064616, 5286.18056848, 5178.46599921, 5073.64292394,
           4996.09317687, 4928.64217002, 4840.85103036, 4777.02248899,
           4734.07319067, 4668.14221256, 4629.6784644 , 4604.04170723,
           4571.57233171, 4521.31548605, 4460.34167673, 4402.72859267,
           4344.38904844, 4292.06954565, 4266.38560573, 4268.05285022,
           4276.05264144, 4271.05051753, 4247.90671997, 4210.51769748,
           4181.10475131, 4168.37295985, 4158.77287475, 4175.66537014,
           4196.54876435, 4205.70600459, 4187.5912329 , 4132.16215453,
           4076.57714652, 4033.44612461, 4003.12294985, 3997.659114  ,
           4005.49455409, 3997.15867066, 3980.80144139, 3952.48902373,
           3913.34122767, 3862.64732013, 3821.88337424, 3829.97002493,
           3841.61061086, 3813.5591777 , 3780.11020795, 3758.81237476,
           3726.80315424, 3701.63067891, 3686.59165315, 3676.71542719,
           3669.57177936, 3681.45398504, 3694.84506278, 3680.13821602,
           3670.78981952, 3673.17354292, 3659.38466996, 3631.02078738,
           3633.40281518, 3628.72992189, 3589.7076796 , 3550.21316276,
           3525.58762807, 3495.28394371, 3478.77031759, 3498.33440902,
    ...
           3395.6703709 , 3415.02232535, 3418.22115408, 3431.49370034,
           3425.29995127, 3422.42361443, 3414.04660616, 3403.52559503,
           3396.53084763, 3387.73173195, 3379.12629815, 3373.38474319,
           3368.55185801, 3364.20044409, 3360.59398   , 3357.64902036,
           3355.09285256, 3353.93538731, 3356.47294464, 3360.37523379,
           3356.86800034, 3341.71626479, 3325.03030682, 3317.00541302,
           3312.00824581, 3299.51168769, 3282.16075266, 3267.2439886 ,
           3257.46847502, 3251.07899002, 3244.64913997, 3239.08900448,
           3234.07062534, 3227.70680978, 3221.27410175, 3214.79795146,
           3209.00781414, 3205.79052014, 3205.75369839, 3207.12839079,
           3204.59770874, 3198.71960138, 3191.07876062, 3182.38339337,
           3174.67817133, 3167.42442929, 3159.30884536, 3151.59404511,
           3147.94278595, 3150.36958969, 3155.35499178, 3154.28407399,
           3143.40365889, 3129.0661275 , 3118.65805078, 3112.63258967,
           3106.30132836, 3099.02010849, 3092.35968023, 3086.83573424,
           3082.68191238, 3080.29880375, 3079.9235029 , 3079.43926301,
           3075.27276504, 3069.42936718, 3068.52037495, 3071.65714631,
           3068.17614263, 3055.97110076, 3045.13693956, 3039.46151169,
           3035.60589776, 3030.59836492, 3025.37443607, 3021.0475862 ,
           3016.18948108, 3010.36748166, 3006.57562897])

  • volume_m3
    (time, dis_along_flowline)
    float64
    1.28e+07 2.056e+07 ... 0.0 0.0

    array([[1.27993428e+07, 2.05559160e+07, 1.14756531e+08, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [1.12936264e+07, 2.15122275e+07, 1.04338325e+08, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [1.03200089e+07, 2.27907901e+07, 9.85235758e+07, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           ...,
           [3.65561670e+06, 2.82687849e+07, 8.32338872e+07, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [3.33761862e+06, 2.79735711e+07, 8.24422603e+07, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [3.50465723e+06, 2.82119124e+07, 8.42577658e+07, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])

  • volume_bsl_m3
    (time, dis_along_flowline)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    description :

    Section volume below sea level

    unit :

    m 3

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]])

  • volume_bwl_m3
    (time, dis_along_flowline)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    description :

    Section volume below water level

    unit :

    m 3

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]])

  • area_m2
    (time, dis_along_flowline)
    float64
    5.402e+05 6.23e+05 ... 0.0 0.0

    array([[ 540202.03506994,  623009.21167763, 3011891.10349759, ...,
                  0.        ,       0.        ,       0.        ],
           [ 537967.55717136,  624235.99407206, 3009122.61128457, ...,
                  0.        ,       0.        ,       0.        ],
           [ 536517.75783621,  625872.41228734, 3007576.31502142, ...,
                  0.        ,       0.        ,       0.        ],
           ...,
           [ 526486.73011012,  632835.73545614, 3003506.58215265, ...,
                  0.        ,       0.        ,       0.        ],
           [ 526003.30801868,  632462.43060294, 3003295.72069481, ...,
                  0.        ,       0.        ,       0.        ],
           [ 526257.29621842,  632763.83598386, 3003779.285302  , ...,
                  0.        ,       0.        ,       0.        ]])

  • thickness_m
    (time, dis_along_flowline)
    float64
    24.12 33.72 38.3 ... 0.0 0.0 0.0

    description :

    Section thickness

    unit :

    m

    array([[24.12456791, 33.7247996 , 38.29592718, ...,  0.        ,
             0.        ,  0.        ],
           [21.33147054, 35.25827759, 34.83531192, ...,  0.        ,
             0.        ,  0.        ],
           [19.51922137, 37.30380036, 32.90244159, ...,  0.        ,
             0.        ,  0.        ],
           ...,
           [ 6.98043671, 46.00795432, 27.8152755 , ...,  0.        ,
             0.        ,  0.        ],
           [ 6.3761591 , 45.54132326, 27.55169868, ...,  0.        ,
             0.        ,  0.        ],
           [ 6.69364435, 45.91807998, 28.15615444, ...,  0.        ,
             0.        ,  0.        ]])

  • ice_velocity_myr
    (time, dis_along_flowline)
    float64
    nan nan nan nan ... 0.0 0.0 0.0 0.0

    description :

    Ice velocity at the surface

    unit :

    m yr-1

    array([[        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [26.4385274 , 31.27999498, 22.11057608, ...,  0.        ,
             0.        ,  0.        ],
           [26.14868126, 31.35046258, 20.81527781, ...,  0.        ,
             0.        ,  0.        ],
           ...,
           [16.91721784, 24.81113009, 18.7257939 , ...,  0.        ,
             0.        ,  0.        ],
           [15.57434689, 23.13726913, 17.93098451, ...,  0.        ,
             0.        ,  0.        ],
           [16.41690876, 24.36672215, 18.97760041, ...,  0.        ,
             0.        ,  0.        ]])

  • calving_bucket_m3
    (time)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    description :

    Flowline calving bucket (volume not yet calved)

    unit :

    m 3

    array([0., 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 nan nan ... 0.0 0.0 0.0 0.0

    description :

    Ice-flux divergence

    unit :

    m yr-1

    array([[        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [-3.40027991,  0.92629545, -4.06779781, ...,  0.        ,
             0.        ,  0.        ],
           [-3.16518522,  0.69258672, -3.28580638, ...,  0.        ,
             0.        ,  0.        ],
           ...,
           [-1.91504792, -1.6554482 , -1.507849  , ...,  0.        ,
             0.        ,  0.        ],
           [-1.84287375, -1.7052272 , -1.50217296, ...,  0.        ,
             0.        ,  0.        ],
           [-1.80853433, -1.74926285, -1.52156382, ...,  0.        ,
             0.        ,  0.        ]])

  • climatic_mb_myr
    (time, dis_along_flowline)
    float64
    nan nan nan nan ... 0.0 0.0 0.0 0.0

    description :

    Climatic mass balance forcing

    unit :

    m yr-1

    array([[       nan,        nan,        nan, ...,        nan,        nan,
                   nan],
           [0.60718254, 0.60718254, 0.60718254, ..., 0.        , 0.        ,
            0.        ],
           [1.35293605, 1.35293605, 1.35293605, ..., 0.        , 0.        ,
            0.        ],
           ...,
           [1.92709701, 1.92709701, 1.92709701, ..., 0.        , 0.        ,
            0.        ],
           [1.23859614, 1.23859614, 1.23859614, ..., 0.        , 0.        ,
            0.        ],
           [2.12601958, 2.12601958, 2.12601958, ..., 0.        , 0.        ,
            0.        ]])

  • dhdt_myr
    (time, dis_along_flowline)
    float64
    nan nan nan nan ... 0.0 0.0 0.0 0.0

    description :

    Thickness change

    unit :

    m yr-1

    array([[        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [-2.79309737,  1.53347799, -3.46061527, ...,  0.        ,
             0.        ,  0.        ],
           [-1.81224917,  2.04552277, -1.93287033, ...,  0.        ,
             0.        ,  0.        ],
           ...,
           [ 0.0120491 ,  0.27164881,  0.41924801, ...,  0.        ,
             0.        ,  0.        ],
           [-0.60427761, -0.46663107, -0.26357682, ...,  0.        ,
             0.        ,  0.        ],
           [ 0.31748525,  0.37675673,  0.60445576, ...,  0.        ,
             0.        ,  0.        ]])

• Indexes: (2)
  • dis_along_flowline
    PandasIndex

    PandasIndex(Index([    0.0,   400.0,   800.0,  1200.0,  1600.0,  2000.0,  2400.0,  2800.0,
            3200.0,  3600.0,
           ...
           67600.0, 68000.0, 68400.0, 68800.0, 69200.0, 69600.0, 70000.0, 70400.0,
           70800.0, 71200.0],
          dtype='float64', name='dis_along_flowline', length=179))

  • time
    PandasIndex

    PandasIndex(Index([2002.0, 2003.0, 2004.0, 2005.0, 2006.0, 2007.0, 2008.0, 2009.0, 2010.0,
           2011.0, 2012.0, 2013.0, 2014.0, 2015.0, 2016.0, 2017.0, 2018.0, 2019.0,
           2020.0],
          dtype='float64', name='time'))

• Attributes: (12)

  calendar :

  365-day no leap

  water_level :

  0

  dx :

  2.0

  creation_date :

  2024-08-25 21:16:01

  fs :

  0

  mb_model_class :

  MultipleFlowlineMassBalance

  map_dx :

  200.0

  class :

  MixedBedFlowline

  glen_a :

  6.42968271645714e-24

  oggm_version :

  1.6.3.dev3+g865475d

  mb_model_hemisphere :

  nh

  description :

  OGGM model output

ds_fl.sel(time=[2003, 2010, 2020]).ice_velocity_myr.plot(hue='time');

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

What’s next?

• return to the OGGM documentation

• back to the table of contents