Creating DTC Glaciers EO-Native Data Cubes (L1)

Creating DTC Glaciers EO-Native Data Cubes (L1)#

In this notebook, we focus on the creation of EO-native data cubes, referred to in the DTC Glaciers project as Level 1 (L1) data cubes.

L1 data cubes form the observational foundation of the DTC Glaciers workflow.
They consist exclusively of Earth Observation (EO) data and in-situ measurements, without the inclusion of model-derived or assimilated quantities. In line with the project design outlined in the proposal, L1 data cubes are generated for individual glaciers, providing a harmonised, analysis-ready representation of the available observations.

All variables within an L1 data cube are projected onto a common spatial grid, temporally aligned where applicable, and pre-processed to enable consistent use in subsequent stages of the DTC, including data assimilation, modelling, and validation workflows (L2 and beyond).

To illustrate the construction and content of L1 data cubes using different EO data sources, this notebook focuses on two contrasting case studies: the Brúarjökull outlet glacier of the Vatnajökull ice cap in Iceland, and Vernagtferner in the Ötztal Alps, Austria.

If required, install the DTCG API using the following command:

!pip install 'dtcg[jupyter] @ git+https://github.com/DTC-Glaciers/dtcg'

Run this command in a notebook cell.

# Imports
import dtcg.integration.oggm_bindings as oggm_bindings
import matplotlib.pyplot as plt

Downloading salem-sample-data...

rgi_id_ice = "RGI60-06.00377"  # Brúarjökull
rgi_id_aut = "RGI60-11.00719"  # Vernagtferner

Enhancing OGGM data cubes with new observations#

CryoTEMPO-EOLIS derived from CryoSat-2#

We developed a method to add four CryoTEMPO-EOLIS variables derived from CryoSat-2 observations to the L1 data cubes.

Variable	Dims	Description
`eolis_gridded_elevation_change`	`(t,y,x)`	Time series of spatial maps of elevation change since January 2011.
`eolis_gridded_elevation_change_sigma`	`(t,y,x)`	Uncertainty (σ) for the gridded elevation-change maps.
`eolis_elevation_change_timeseries`	`(t)`	Spatially aggregated 1D elevation-change series since January 2011.
`eolis_elevation_change_sigma_timeseries`	`(t)`	Uncertainty (σ) for the aggregated series.

To use this on your own laptop, you need to register on the Specklia website and obtain an API key. You then need to provide this key to DTC Glaciers in the following way:

# import os
# os.environ["SPECKLIA_API_KEY"] = ENTER YOUR SPECKLIA API KEY HERE

You can add the data using add_data together with DatacubeCryotempoEolis. This step can take some time, as up-to-date data are downloaded from the CryoTEMPO-EOLIS server and reprojected to the L1 data cubes.

Let’s visualise these new data:

../_images/2537e203688eca1d20b9f25ad0a0cae09727fa8df2d093f5fe0ec160c58b02ce.png

../_images/eb5e71d20f9cd83520c570efefdbf3e256a385090fc577a713c37059a21202bf.png

Glacier surface classification from Sentinel-2#

For DTC Glaciers, ENVEO provides a 10 m gridded glacier surface classification product derived from Sentinel-2 data. This dataset can be added to the L1 data cubes in the same way as shown above.

Variable	Dims	Description
`sfc_type_data`	`(t_sfc_type,y,x)`	Glacier surface classification.
`sfc_type_uncertainty`	`(t_sfc_type,y,x)`	Glacier surface classification uncertainty.
`sfc_type_snowline`	`(snowcover_frac, t_sfc_type)`	Snowline altitude derived from glacier surface classification for three different snow cover area fractions (25%, 50% and 75%) per 30m elevation band.

The snowline was derived within DTC Glaciers for model validation and could also be assimilated in the future. To derive it, we aggregated the gridded glacier surface classification into 30 m elevation bands and calculated the snow cover area fraction for each band, excluding cloud-covered pixels.

We then identified the lowest elevation band where one of the three thresholds - 25%, 50%, or 75% - was exceeded for the first time.

Let’s look at one example:

Show code cell source

Hide code cell source

import numpy as np
import xarray as xr
from dtcg.datacube.surface_type import (
    get_categories_per_elevation_band,
    exclude_empty_elevation_bins,
    exclude_dates_with_too_much_cloud_cover,
    get_snowline,
)


def plot_surface_classification_per_elev_band(datacube_l1, date):
    ds_obs = get_categories_per_elevation_band(datacube_l1)
    ds_obs = exclude_empty_elevation_bins(ds=ds_obs)
    ds_obs = exclude_dates_with_too_much_cloud_cover(ds=ds_obs)

    ds_use = ds_obs.sel(t_sfc_type=date)
    ds_use = xr.decode_cf(ds_use)

    bin_idx, bin_elev = get_snowline(ds_use)

    fig, axs = plt.subplots(1, 3, figsize=(13, 5), gridspec_kw={"wspace": 0.1})

    ax_total = axs[0]
    ax_relative = axs[1]
    ax_uncertainty = axs[2]

    def plot_horizontal_stacked_bar(ax, ds, categories, var, normalize=False):
        category_dict = eval(
            ds[categories].attrs["code"]
        )  # ast.literal_eval(ds[categories].attrs['code'])
        previous_values = None
        for category in ds_use[categories]:
            current_values = ds.loc[{categories: category}][var].values
            if normalize:
                number_grid_points_elev_bin = ds.category_counts.sum(dim="category")
                current_values = (
                    current_values.astype(float) / number_grid_points_elev_bin
                )
            ax.barh(
                ds.elevation_bin,
                current_values,
                left=previous_values,
                label=category_dict[int(category)],
            )

            if previous_values is None:
                previous_values = current_values
            else:
                previous_values = previous_values + current_values

        ax.invert_yaxis()

    plot_horizontal_stacked_bar(ax_total, ds_use, "category", "category_counts")
    plot_horizontal_stacked_bar(
        ax_relative, ds_use, "category", "category_counts", normalize=True
    )
    plot_horizontal_stacked_bar(
        ax_uncertainty, ds_use, "uncertainty_flag", "uncertainty_counts", normalize=True
    )

    def add_snowline(ax):
        if np.isneginf(bin_elev[1]):
            extra_label = "fully snow covered"
        elif np.isposinf(bin_elev[1]):
            extra_label = "fully snow free"
        else:
            extra_label = f"{bin_elev[1]:.0f} m"
        ax.axhline(bin_idx[0], color="k", linestyle="--")
        ax.axhline(
            bin_idx[1], color="k", linestyle="-", label=f"snowline ({extra_label})"
        )
        ax.axhline(bin_idx[2], color="k", linestyle="--", label="snowline uncertainty")

    def set_legend(ax):
        ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.05))

    ax_total.set_yticks(
        ds_use.elevation_bin, [int(elev) for elev in ds_use.lower_elevation.values]
    )
    ax_total.set_ylabel("Altitude [m]")
    ax_total.set_title("category count")
    add_snowline(ax_total)
    set_legend(ax_total)

    ax_relative.set_title("category count relative")
    ax_relative.set_yticklabels([])
    add_snowline(ax_relative)
    set_legend(ax_relative)

    ax_uncertainty.set_title("uncertainty count relative")
    ax_uncertainty.set_yticklabels([])
    add_snowline(ax_uncertainty)
    set_legend(ax_uncertainty)

    fig.suptitle(
        f"{datacube_l1.attrs['RGI-ID']}, {np.datetime_as_string(ds_use.t_sfc_type.values, unit='D')}"
    )

    plt.show()

plot_surface_classification_per_elev_band(
    datacube_l1=dtcg_oggm_aut.datacube_l1, date=dtcg_oggm_aut.datacube_l1.t_sfc_type[1]
)

../_images/dbbf48ba2e1afe001cb00ee8540d54ae1a33ca8c0072128388ce57080981c6ba.png

Let’s visualize the derived snowline on top of the original gridded data:

Show code cell source

Hide code cell source

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import seaborn as sns
import matplotlib.lines as mlines


def build_categorical_colormap(class_labels, palette_name="colorblind"):

    # Extract integer classes and sort them
    classes = sorted(
        [k for k in class_labels.keys() if isinstance(k, (int, np.integer))]
    )
    n = len(classes)

    # Pick n colors from seaborn palette
    colors = sns.color_palette(palette_name, n)

    # Create colormap
    cmap = mcolors.ListedColormap(colors)

    # Set NaN color ("bad" values)
    # cmap.set_bad("lightgray")

    # Boundaries between classes — last boundary is +1
    bounds = [c - 0.5 for c in classes] + [classes[-1] + 0.5]
    norm = mcolors.BoundaryNorm(bounds, cmap.N)

    return norm, cmap, classes


def plot_gridded_surface_classification(
    datacube_l1, date, x_zoom=70, y_zoom=50, snowline=None
):

    ds_plot = datacube_l1.isel(
        x=slice(x_zoom, len(datacube_l1.x) - x_zoom),
        y=slice(y_zoom, len(datacube_l1.y) - y_zoom),
    )

    ds_plot = ds_plot.sel(t_sfc_type=date)
    ds_plot = xr.decode_cf(ds_plot)

    # Integer classes in your dataset
    dict_labels_1 = {
        0: "unclassified",
        1: "snow",
        2: "firn \n old snow \n bright ice",
        3: "clean ice",
        4: "debris",
        5: "cloud",
        "nan": "no data",
    }  # ast.literal_eval(ds_plot.sfc_type_data.code)

    dict_labels_2 = {
        1: "low uncertainty \n illuminated pixel",
        2: "medium uncertainty \n illuminated pixel",
        3: "high uncertainty \n illuminated pixel",
        5: "cloud",
        11: "low uncertainty \n shaded pixel",
        12: "medium uncertainty \n shaded pixel",
        13: "high uncertainty \n shaded pixel",
        "nan": "no data",
    }  # ast.literal_eval(ds_plot.sfc_type_uncertainty.code)

    norm_1, cmap_1, classes_1 = build_categorical_colormap(dict_labels_1)
    norm_2, cmap_2, classes_2 = build_categorical_colormap(dict_labels_2)

    # Your data array
    da_1 = ds_plot.where(ds_plot.glacier_mask).sfc_type_data
    da_2 = ds_plot.where(ds_plot.glacier_mask).sfc_type_uncertainty

    fig, axs = plt.subplots(1, 2, figsize=(10, 5))

    ax_1 = axs[0]
    ax_2 = axs[1]

    mappable_1 = da_1.plot(ax=ax_1, cmap=cmap_1, norm=norm_1, add_colorbar=False)
    mappable_2 = da_2.plot(ax=ax_2, cmap=cmap_2, norm=norm_2, add_colorbar=False)

    # add derived snowline

    cs = ds_plot.where(ds_plot.glacier_mask).topo_smoothed.plot.contour(
        ax=ax_1, levels=[snowline], colors="k", linewidths=0.5
    )
    proxy = mlines.Line2D(
        [], [], color="k", linewidth=0.5, label=f"derived snowline at {snowline} m"
    )
    ax_1.legend(handles=[proxy])

    cs = ds_plot.where(ds_plot.glacier_mask).topo_smoothed.plot.contour(
        ax=ax_2, levels=[snowline], colors="k", linewidths=0.5
    )
    proxy_2 = mlines.Line2D(
        [], [], color="k", linewidth=0.5, label=f"derived snowline at {snowline} m"
    )
    ax_2.legend(handles=[proxy_2])

    # Add labeled colorbar
    cbar_1 = plt.colorbar(
        mappable_1, ax=ax_1, orientation="horizontal", ticks=classes_1
    )
    cbar_1.ax.set_xticklabels([f"{dict_labels_1[c]}" for c in classes_1], rotation=90)
    cbar_1.set_label(ds_plot.sfc_type_data.long_name)

    cbar_2 = plt.colorbar(
        mappable_2,
        ax=ax_2,
        ticks=[1, 2, 3.5, 7.5, 11, 12, 13],
        orientation="horizontal",
    )  # classes_2)
    cbar_2.ax.set_xticklabels([f"{dict_labels_2[c]}" for c in classes_2], rotation=90)
    cbar_2.set_label(ds_plot.sfc_type_uncertainty.long_name)

    ax_1.axis("equal")
    ax_2.axis("equal")

    ax_1.set_title(
        f"{datacube_l1.attrs['RGI-ID']}, {np.datetime_as_string(da_1.t_sfc_type.values, unit='D')}"
    )
    ax_2.set_title(
        f"{datacube_l1.attrs['RGI-ID']}, {np.datetime_as_string(da_2.t_sfc_type.values, unit='D')}"
    )

    fig.tight_layout()
    plt.show()

plot_gridded_surface_classification(
    datacube_l1=dtcg_oggm_aut.datacube_l1,
    date=dtcg_oggm_aut.datacube_l1.t_sfc_type[1],
    snowline=3120,
)

../_images/e273deb675da4b8ee11f07b3b4a9a37bd91951efa32306cc4bbd5934fb54ec98.png

In-situ mass balance observations from WGMS#

For validation, the DTCG API can add in-situ mass balance observations from WGMS to the L1 data cubes for glaciers where such observations are available, using the same approach as shown above.

Variable	Dims	Description
`wgms_mb`	`(t_wgms)`	Specific mass balance observation as reported tothe World Glacier Monitoring Service.
`wgms_mb_unc`	`(t_wgms)`	Specific mass balance observation uncertainty as reported to the World Glacier Monitoring Service. If no value was reported it is set to 200 mm w.e..

plt.errorbar(
    x=dtcg_oggm_ice.datacube_l1.t_wgms,
    y=dtcg_oggm_ice.datacube_l1.wgms_mb,
    yerr=dtcg_oggm_ice.datacube_l1.wgms_mb_unc,
    fmt=".--",
    capsize=2,
)
plt.grid("on")
plt.ylabel("Specific MB (mm w.e. yr-1)")
plt.xlabel("Year")
plt.title(dtcg_oggm_ice.datacube_l1.attrs["RGI-ID"])
plt.show()

../_images/bd73ab1eaf3712c2be483020e24c1cdc44639c140ac40ce6fdcdb2f012f79733.png

Preprocessed data cubes for case study regions#

Preprocessed L1 data cubes, including all currently supported observations, are available for the case study regions in Iceland and Austria. You can obtain the URL of the latest version via DEFAULT_L1_DATACUBE_URL and use the DTCG API’s dedicated data cube streaming functionality as shown below:

from dtcg import DEFAULT_L1_DATACUBE_URL
from dtcg.api.external.call import StreamDatacube

streamer = StreamDatacube(server=DEFAULT_L1_DATACUBE_URL)
datacube = streamer.stream_datacube(glacier=rgi_id_ice)
datacube

<xarray.DataTree>
Group: /
└── Group: /L1
        Dimensions:                                  (y: 370, x: 436, t: 177, t_wgms: 32)
        Coordinates:
          * y                                        (y) float32 1kB 7.201e+06 ... 7....
          * x                                        (x) float32 2kB 3.94e+05 ... 4.8...
          * t                                        (t) datetime64[ns] 1kB 2011-01-1...
          * t_wgms                                   (t_wgms) int64 256B 1993 ... 2024
        Data variables: (12/21)
            consensus_ice_thickness                  (y, x) float32 645kB dask.array<chunksize=(370, 436), meta=np.ndarray>
            eolis_elevation_change_sigma_timeseries  (t) float64 1kB dask.array<chunksize=(177,), meta=np.ndarray>
            eolis_elevation_change_timeseries        (t) float64 1kB dask.array<chunksize=(177,), meta=np.ndarray>
            eolis_gridded_elevation_change           (t, y, x) float64 228MB dask.array<chunksize=(177, 60, 60), meta=np.ndarray>
            eolis_gridded_elevation_change_sigma     (t, y, x) float64 228MB dask.array<chunksize=(177, 60, 60), meta=np.ndarray>
            glacier_ext                              (y, x) int8 161kB dask.array<chunksize=(370, 436), meta=np.ndarray>
            ...                                       ...
            spatial_ref                              int64 8B ...
            topo                                     (y, x) float32 645kB dask.array<chunksize=(370, 436), meta=np.ndarray>
            topo_smoothed                            (y, x) float32 645kB dask.array<chunksize=(370, 436), meta=np.ndarray>
            topo_valid_mask                          (y, x) int8 161kB dask.array<chunksize=(370, 436), meta=np.ndarray>
            wgms_mb                                  (t_wgms) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
            wgms_mb_unc                              (t_wgms) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
        Attributes:
            Conventions:         CF-1.12
            comment:             The DTC Glaciers project is developed under the Euro...
            date_created:        2026-01-10T23:48:41.290684
            RGI-ID:              RGI60-06.00377
            glacier_attributes:  {'rgi_id': 'RGI60-06.00377', 'glims_id': 'G343733E64...
            title:               Datacube of glacier-domain variables.
            summary:             Resampled glacier-domain variables from multiple sou...

Or you can use xr.open_datatree for data streaming:

def get_l1_data_tree(rgi_id):
    return xr.open_datatree(
        f"{DEFAULT_L1_DATACUBE_URL}{rgi_id}.zarr",
        chunks={},
        engine="zarr",
        consolidated=True,
        decode_cf=True,
    )

Streaming means that only a view of the available data, together with the corresponding coordinates, is provided first. The actual data are only downloaded when needed. This approach saves bandwidth and makes it fast to explore large datasets.

The preprocessed data cubes (L1 and L2, explained in this notebook) are currently available for all RGI6 glaciers of Vatnajökull (view in OpenStreetMap):

../_images/7f9376123a29ecd7064006d1a337793f9bf864e37774284257a2db98ba472efe.png

And for the Ötztal region in the Austrian Alps (view in OpenStreetMap):

../_images/c1f69575c210a4829311198fa2d9f7de27c1cb62d8763bfc012a557f9a4b458f.png

To obtain the RGI6 ID of a glacier, you can use the GLIMS Viewer. Enable the RGIv6 IDs in the menu in the top-right corner, zoom to the glacier of interest, and click on its outline to display the metadata.