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.

NOTE: This code needs OGGM v1.5.3 to run.

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

2023-03-07 12:28:46: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-07 12:28:46: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-07 12:28:46: oggm.cfg: Multiprocessing: using all available processors (N=2)

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.4/L3-L5_files/CRU/elev_bands/qc3/pcp2.5/no_match/'

gdirs = workflow.init_glacier_directories(rgi_ids,
                                          from_prepro_level=from_prepro_level,
                                          prepro_base_url=base_url,
                                          prepro_border=prepro_border)

2023-03-07 12:28:46: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2023-03-07 12:28:46: 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>
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.299e+03
    topo_smoothed    (y, x) float32 5.113e+03 5.019e+03 ... 5.291e+03 5.295e+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

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.299e+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.],
           ...,
           [3974., 3866., 3783., ..., 5252., 5260., 5263.],
           [4025., 3917., 3812., ..., 5279., 5286., 5287.],
           [4097., 4005., 3916., ..., 5291., 5291., 5299.]], dtype=float32)

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

    units :

    m

    long_name :

    DEM topography smoothed with radius: 3e+02 m

    array([[5113.0107, 5019.3696, 4943.0776, ..., 5468.2134, 5543.3647,
            5599.1377],
           [5138.9595, 5046.7725, 4977.5405, ..., 5396.9126, 5451.853 ,
            5503.463 ],
           [5128.379 , 5063.5703, 5022.3716, ..., 5261.8374, 5320.353 ,
            5378.3765],
           ...,
           [3943.0627, 3864.9866, 3780.063 , ..., 5257.0786, 5258.0728,
            5262.429 ],
           [4007.28  , 3926.4658, 3832.7732, ..., 5287.942 , 5280.9517,
            5283.8403],
           [4061.5234, 3986.9854, 3895.4531, ..., 5302.535 , 5291.0063,
            5295.076 ]], 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(Float64Index([-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='float64', name='x', length=479))

  • y
    PandasIndex

    PandasIndex(Float64Index([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='float64', 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();

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

2023-03-07 12:28:59: 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: ITS-live

This is an example on how to extract velocity fields from the ITS_live Regional Glacier and Ice Sheet Surface Velocities Mosaic (Gardner, A. et al 2019) at 120 m resolution and reproject this data to the OGGM-glacier grid. This only works where ITS-live data is available! (not in the Alps).

The data source used is https://its-live.jpl.nasa.gov/#data.

Currently the only data downloaded is the 120m composite for both (u, v) and their uncertainty. The composite is computed from the 1985 to 2018 average. If you want more velocity products, feel free to open a new topic on the OGGM issue tracker!

this will download severals large dataset (2~800MB) depending on your connection might take a bit …

from oggm.shop import its_live
workflow.execute_entity_task(its_live.velocity_to_gdir, gdirs);

2023-03-07 12:29:11: oggm.workflow: Execute entity tasks [velocity_to_gdir] on 1 glaciers

By applying the entity task its_live.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.

According to the ITS_LIVE documentation velocities are given in ground units (i.e. absolute velocities). We then use bilinear interpolation to reproject the velocities to the local glacier map by re-projecting the vector distances (For more details feel free to checl out the code).

By specifying add_error=True on its_live.velocity_to_gdir(gdir, add_error=True), we also reproject and scale the error for each component (evx, evy).

Now we can read in all the gridded data that comes with OGGM, including the ITS_Live 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.299e+03
    topo_smoothed            (y, x) float32 5.113e+03 5.019e+03 ... 5.295e+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
    obs_icevel_x             (y, x) float32 -0.3372 -0.3311 ... -16.64 -19.13
    obs_icevel_y             (y, x) float32 1.797 5.185 8.147 ... 33.85 38.78
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: (8)
  • topo
    (y, x)
    float32
    5.128e+03 4.999e+03 ... 5.299e+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.],
           ...,
           [3974., 3866., 3783., ..., 5252., 5260., 5263.],
           [4025., 3917., 3812., ..., 5279., 5286., 5287.],
           [4097., 4005., 3916., ..., 5291., 5291., 5299.]], dtype=float32)

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

    units :

    m

    long_name :

    DEM topography smoothed with radius: 3e+02 m

    array([[5113.0107, 5019.3696, 4943.0776, ..., 5468.2134, 5543.3647,
            5599.1377],
           [5138.9595, 5046.7725, 4977.5405, ..., 5396.9126, 5451.853 ,
            5503.463 ],
           [5128.379 , 5063.5703, 5022.3716, ..., 5261.8374, 5320.353 ,
            5378.3765],
           ...,
           [3943.0627, 3864.9866, 3780.063 , ..., 5257.0786, 5258.0728,
            5262.429 ],
           [4007.28  , 3926.4658, 3832.7732, ..., 5287.942 , 5280.9517,
            5283.8403],
           [4061.5234, 3986.9854, 3895.4531, ..., 5302.535 , 5291.0063,
            5295.076 ]], 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)

  • obs_icevel_x
    (y, x)
    float32
    -0.3372 -0.3311 ... -16.64 -19.13

    units :

    m yr-1

    long_name :

    ITS LIVE velocity data in x map direction

    array([[-3.3717626e-01, -3.3110973e-01, -3.3433609e+00, ...,
            -1.0667664e-02, -1.0126008e-01, -3.4660590e-01],
           [ 5.3639662e-02,  1.2223738e+00,  1.2713091e-01, ...,
             9.7065814e-02,  2.1466576e-02, -2.1978757e-01],
           [ 3.7502187e-01,  1.2963842e+00,  9.0304291e-01, ...,
             2.0853087e-01,  8.7231494e-02,  7.5737543e-02],
           ...,
           [-4.1035689e-02, -7.6566137e-02, -7.3401690e-02, ...,
            -1.3299100e+01, -1.7874100e+01, -1.6857393e+01],
           [-1.0310390e-02, -1.9115219e-01, -1.7953111e-01, ...,
            -1.4452621e+01, -2.0741724e+01, -1.9713362e+01],
           [-8.3360513e-04, -1.0499098e-01, -1.6805069e-01, ...,
            -1.1871961e+01, -1.6640432e+01, -1.9127012e+01]], dtype=float32)

  • obs_icevel_y
    (y, x)
    float32
    1.797 5.185 8.147 ... 33.85 38.78

    units :

    m yr-1

    long_name :

    ITS LIVE velocity data in y map direction

    array([[ 1.79682636e+00,  5.18474722e+00,  8.14693832e+00, ...,
            -4.86147553e-01, -3.16282541e-01, -4.03148592e-01],
           [ 2.07681346e+00,  5.12276125e+00,  6.94922876e+00, ...,
            -2.11663038e-01, -7.33599141e-02, -4.68726218e-01],
           [ 2.75870657e+00,  4.72603559e+00,  6.15555000e+00, ...,
             1.12408385e-01,  3.07224423e-01,  3.72547358e-01],
           ...,
           [ 1.57866269e-01,  1.26633152e-01,  6.88849315e-02, ...,
             3.74994545e+01,  4.64094391e+01,  4.40526390e+01],
           [ 2.72619762e-02,  1.95747986e-02, -3.38266231e-02, ...,
             2.98384285e+01,  4.55910683e+01,  4.56010513e+01],
           [-7.31839566e-03, -8.27954710e-02, -9.56750214e-02, ...,
             1.77792721e+01,  3.38503227e+01,  3.87761536e+01]], dtype=float32)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Float64Index([-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='float64', name='x', length=479))

  • y
    PandasIndex

    PandasIndex(Float64Index([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='float64', 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);

# get the velocity data
u = ds.obs_icevel_x.where(ds.glacier_mask == 1)
v = ds.obs_icevel_y.where(ds.glacier_mask == 1)
ws = (u**2 + v**2)**0.5

The ds.glacier_mask == 1 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, 100, '100 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', 'obs_icevel_x', 'obs_icevel_y'],
                              preserve_totals=[True, False, False]  # I"m actually not sure if preserving totals is meaningful with velocities - likely not
                              )

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	obs_icevel_x	obs_icevel_y	bin_elevation	dx	width
dis_along_flowline
60.172585	40000.0	8345.080078	0.244304	18.147763	-0.191088	0.212503	8355.0	120.345169	332.377280
140.404661	80000.0	8235.238281	0.642076	13.891147	-0.119776	0.298178	8235.0	40.118984	1994.068447
173.940860	40000.0	8197.490234	0.838840	21.082249	-0.173281	0.277646	8205.0	26.953414	1484.042070
202.717642	40000.0	8150.462891	0.775495	25.654437	0.060681	0.333409	8145.0	30.600150	1307.183150
236.363889	40000.0	8075.439941	0.685386	14.702222	0.089171	0.122661	8085.0	36.692346	1090.145619
...	...	...	...	...	...	...	...	...	...
42308.885990	760000.0	3525.279541	0.104131	139.481158	-4.474717	-2.795007	3525.0	287.057248	2647.555512
42583.856935	400000.0	3495.227295	0.113627	101.916241	-3.964900	-2.077398	3495.0	262.884641	1521.579952
42814.941987	200000.0	3472.789551	0.149416	98.673038	-4.508717	-1.815381	3465.0	199.285463	1003.585494
43010.329233	200000.0	3441.432373	0.155404	59.133279	-2.572173	-1.226399	3435.0	191.489029	1044.446259
43224.984756	120000.0	3407.177490	0.125482	32.996786	-0.836834	-0.399884	3405.0	237.822018	504.579018

157 rows × 9 columns

df['obs_velocity'] = (df['obs_icevel_x']**2 + df['obs_icevel_y']**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', 'obs_icevel_x', 'obs_icevel_y'],
                                       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	obs_icevel_x	obs_icevel_y
dis_along_flowline
200.0	1350.508131	5.402033e+05	24.852592	0.038586	0.328143
600.0	1557.488984	6.229956e+05	12.058255	-0.115721	0.288789
1000.0	7529.527187	3.011811e+06	32.532095	2.174206	0.080531
1400.0	15906.311989	6.362525e+06	29.213418	0.595236	0.310646
1800.0	35987.148540	1.439486e+07	42.631021	0.415668	0.740916
...	...	...	...	...	...
41400.0	2128.789125	8.515156e+05	188.005470	-5.837341	-11.296937
41800.0	1795.906259	7.183625e+05	163.145085	-7.081799	-7.600463
42200.0	2387.409586	9.549638e+05	148.023098	-5.238625	-3.783761
42600.0	1515.264571	6.061058e+05	100.197741	-4.002890	-2.059094
43000.0	1063.245975	4.252984e+05	60.325326	-2.674549	-1.257536

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['obs_icevel_x']**2 + df_regular['obs_icevel_y']**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',
 'inversion_glen_a': 2.2535112971797524e-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']

2023-03-07 12:29:30: oggm.cfg: PARAMS['inversion_glen_a'] changed from `2.4e-24` to `2.2535112971797524e-24`.

# Since we overwrote the inversion flowlines we have to do stuff again
utils.apply_test_ref_tstars()
tasks.local_t_star(gdir)
tasks.mu_star_calibration(gdir)
tasks.prepare_for_inversion(gdir, add_debug_var=True)
tasks.mass_conservation_inversion(gdir)
tasks.compute_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();

Velocities during the 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

# if you start to look into velocities, you will see that OGGM is quite large 
# regarding numerical stabilities - here is a better default (yet not perfect)
cfg.PARAMS['cfl_number'] = 0.01

2023-03-07 12:29:30: oggm.cfg: PARAMS['store_fl_diagnostics'] changed from `False` to `True`.
2023-03-07 12:29:30: oggm.cfg: PARAMS['cfl_number'] changed from `0.02` to `0.01`.

tasks.run_from_climate_data(gdir);

/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 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>
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:
    point_lons          (dis_along_flowline) float64 75.88 75.88 ... 76.66 76.67
    point_lats          (dis_along_flowline) float64 36.07 36.07 ... 36.07 36.07
    bed_h               (dis_along_flowline) float64 8.129e+03 ... 3.007e+03
    volume_m3           (time, dis_along_flowline) float64 1.377e+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
    area_m2             (time, dis_along_flowline) float64 5.402e+05 ... 0.0
    thickness_m         (time, dis_along_flowline) float64 25.99 36.33 ... 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
Attributes:
    class:                MixedBedFlowline
    map_dx:               200.0
    dx:                   2.0
    description:          OGGM model output
    oggm_version:         1.5.3
    calendar:             365-day no leap
    creation_date:        2023-03-07 12:29:30
    water_level:          0
    glen_a:               2.2535112971797524e-24
    fs:                   0
    mb_model_class:       MultipleFlowlineMassBalance
    mb_model_hemisphere:  nh

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: (10)
  • 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.129e+03 7.691e+03 ... 3.007e+03

    description :

    Bed elevation along the flowline

    unit :

    m

    array([8128.91513405, 7691.13170378, 7510.58633852, 7246.90520026,
           7042.71036754, 6796.40588114, 6558.3301199 , 6358.46386639,
           6097.94300963, 5878.2440223 , 5683.15841692, 5518.01019407,
           5398.62471699, 5283.84894279, 5176.30753594, 5071.40370433,
           4993.95641361, 4926.78428474, 4838.9477781 , 4775.20168075,
           4732.55346493, 4666.59517229, 4628.31828239, 4602.93430046,
           4570.63769102, 4520.40443824, 4459.39736486, 4401.79851366,
           4343.35104746, 4290.87824532, 4265.3015801 , 4267.1936893 ,
           4275.43340172, 4270.65100458, 4247.51719857, 4210.07289469,
           4180.74644965, 4168.13310843, 4158.65602363, 4175.70442438,
           4196.74410656, 4206.0380911 , 4188.15155546, 4132.8409503 ,
           4077.41298782, 4034.4506279 , 4004.22241208, 3998.78923332,
           4006.72984887, 3998.57875546, 3982.4954052 , 3954.51968769,
           3915.69212487, 3865.3830426 , 3824.94899468, 3833.03499297,
           3844.83514265, 3817.26408869, 3784.37606942, 3763.56907305,
           3732.13911914, 3707.45831852, 3692.79592611, 3683.30927909,
           3676.58870748, 3688.66275124, 3701.82496469, 3687.12997152,
           3677.96963157, 3680.75349881, 3667.65938715, 3639.81781685,
           3642.11863899, 3637.21494281, 3598.47333174, 3559.57927028,
           3535.49247871, 3505.94589459, 3490.10637906, 3509.57832412,
    ...
           3415.27861241, 3432.91886156, 3433.88783344, 3446.05325835,
           3419.12957665, 3403.41930471, 3378.11034317, 3348.08896802,
           3396.53145556, 3387.73208598, 3379.1264285 , 3373.38492634,
           3368.55194186, 3364.2002531 , 3360.59373743, 3357.64888115,
           3355.09291753, 3353.93567926, 3356.47342392, 3360.37573296,
           3356.86844231, 3341.71679679, 3325.0308002 , 3317.00568497,
           3312.00842881, 3299.51186469, 3282.16085483, 3267.24407544,
           3257.46871454, 3251.07948217, 3244.6497637 , 3239.08948298,
           3234.07079021, 3227.70662272, 3221.27355922, 3214.79725169,
           3209.00738083, 3205.79057165, 3205.75402596, 3207.12878622,
           3204.59802108, 3198.71971913, 3191.07880376, 3182.38354826,
           3174.67824984, 3167.42417347, 3159.30851391, 3151.59407422,
           3147.94336266, 3150.37052749, 3155.35601543, 3154.28496026,
           3143.40443906, 3129.06680039, 3118.65858664, 3112.63305909,
           3106.30164701, 3099.02031129, 3092.35997413, 3086.83610516,
           3082.68230023, 3080.29905511, 3079.9235937 , 3079.43954113,
           3075.2731983 , 3069.4297313 , 3068.5208008 , 3071.65772001,
           3068.17681682, 3055.97173942, 3045.13726132, 3039.46150185,
           3035.60584233, 3030.59837538, 3025.37448824, 3021.04768202,
           3016.18962454, 3010.3676768 , 3006.57581382])

  • volume_m3
    (time, dis_along_flowline)
    float64
    1.377e+07 2.211e+07 ... 0.0 0.0

    array([[1.37691489e+07, 2.21075405e+07, 1.23548651e+08, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [1.30237531e+07, 2.26988913e+07, 1.18155884e+08, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [1.24372892e+07, 2.34222549e+07, 1.14117703e+08, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           ...,
           [5.41775817e+06, 2.76052760e+07, 8.10259141e+07, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [5.06659289e+06, 2.76832512e+07, 7.97724704e+07, ...,
            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
           [4.81196395e+06, 2.78402614e+07, 7.90597815e+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([[ 540203.25066458,  622995.67956796, 3011810.98497579, ...,
                  0.        ,       0.        ,       0.        ],
           [ 539098.24595405,  623754.58156288, 3010378.21245773, ...,
                  0.        ,       0.        ,       0.        ],
           [ 538227.25359917,  624681.64679929, 3009304.88518565, ...,
                  0.        ,       0.        ,       0.        ],
           ...,
           [ 527690.55973757,  630015.86768656, 3000494.79742773, ...,
                  0.        ,       0.        ,       0.        ],
           [ 527157.910302  ,  630114.87356622, 3000160.58227769, ...,
                  0.        ,       0.        ,       0.        ],
           [ 526771.35085701,  630314.18367405, 2999970.536073  , ...,
                  0.        ,       0.        ,       0.        ]])

  • thickness_m
    (time, dis_along_flowline)
    float64
    25.99 36.33 41.25 ... 0.0 0.0 0.0

    description :

    Section thickness

    unit :

    m

    array([[25.98895742, 36.33346005, 41.24733899, ...,  0.        ,
             0.        ,  0.        ],
           [24.60770153, 37.28208754, 39.45637334, ...,  0.        ,
             0.        ,  0.        ],
           [23.51896109, 38.44091909, 38.11471425, ...,  0.        ,
             0.        ,  0.        ],
           ...,
           [10.34809376, 45.1086952 , 27.10210456, ...,  0.        ,
             0.        ,  0.        ],
           [ 9.68228196, 45.23245255, 26.68433562, ...,  0.        ,
             0.        ,  0.        ],
           [ 9.19908266, 45.48159018, 26.44677786, ...,  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],
           [13.45460832, 15.80200126, 11.53235953, ...,  0.        ,
             0.        ,  0.        ],
           [13.29587888, 15.71897178, 11.02320889, ...,  0.        ,
             0.        ,  0.        ],
           ...,
           [ 7.38774971,  9.92578638,  6.37266054, ...,  0.        ,
             0.        ,  0.        ],
           [ 7.06527363,  9.58383167,  6.19022663, ...,  0.        ,
             0.        ,  0.        ],
           [ 6.90502193,  9.44300802,  6.12126594, ...,  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.])

• Indexes: (2)
  • dis_along_flowline
    PandasIndex

    PandasIndex(Float64Index([    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(Float64Index([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)

  class :

  MixedBedFlowline

  map_dx :

  200.0

  dx :

  2.0

  description :

  OGGM model output

  oggm_version :

  1.5.3

  calendar :

  365-day no leap

  creation_date :

  2023-03-07 12:29:30

  water_level :

  0

  glen_a :

  2.2535112971797524e-24

  fs :

  0

  mb_model_class :

  MultipleFlowlineMassBalance

  mb_model_hemisphere :

  nh

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