EO-LINCS project: Cube generation for Scientific Case Study (SCS) 2¶

Forest recovery post disturbance¶

Objective: The SCS2 aims to estimate canopy height and above-ground biomass maps using deep learning methods based on remote sensing data. The canopy height and above-ground biomass maps are subsequently used to calculate recovery curves, which in turn can be used to estimate carbon budget.

Outcomes: New high-resolution height/biomass maps that are expected to enable the monitoring of biomass at finer scales, in particular the impact of fine scale forest disturbances due to management practices such as thinning and the impact of natural disturbances (insects attacks, droughts, fires and windthrown in regions of interest). Analysis of forest recovery depending on environmental factors (such as climate, soil composition and pH) and the nature and intensity of disturbance shall aid in the optimization of forest management considering potential increased future disturbances.

Required Datasets:

• Global Age Mapping Integration (GAMI) dataset
• European Forest Disturbance Atlas
• Canopy height and biomass map for Europe
• Copernicus Tree cover density 2015
• Copernicus Dominant Leaf Type 2015
• Copernicus Forest Type 2015

---

The following notebook shows how the users can load data from various sources defined in scs2_config.yml using the MultiSourceDataStore tool.

What You Can Do with This Notebook¶

• Load datasets from various sources as defined in the scs2_config.yml
• View the progress of each data request to the MultiSourceDataStore
• Quickly preview the datasets by downsampling and plotting them.

Requirements¶

Before you begin, ensure that all necessary dependencies are installed and that you have generated API token credentials for CLMS:

• Install xcube-multistore via conda-forge by running: conda install --channel conda-forge xcube-multistore
• Generate CLMS API token credentials by following the offical instructions. Save the credentials in a file named clms-credentials.json, placing it in the same directory as this notebook. Note that the filename can be customized in the scs2_config.yml configuration file.

Once you have it installed, you are ready to proceed.

This Multistore mainly works with a file called scs2_config.yml which is at the same file level as this notebook. To understand what goes into the schema, you can read more here.

---

Let's import the MultiSourceDataStore

In [1]:

from xcube_multistore import MultiSourceDataStore

from xcube_multistore import MultiSourceDataStore

You can find out how to fill out the config file by also using this super helpful function get_config_schema(). Run it and try expand the fields to learn more about the possible properties that the configuration file accepts along with the Configuration Guide.

In [2]:

MultiSourceDataStore.get_config_schema()

MultiSourceDataStore.get_config_schema()

Out[2]:

<xcube.util.jsonschema.JsonObjectSchema at 0x7e1b612951d0>

To set up the config file, please refer to the notebook setup_config.ipynb or the Configuration Guide in the documentation.

We can view the configuration using the method display_config, executing the follwoing cell. It shows, that for the user-defined dataset ID clms and senf, some data fusion of multiple data sources will be performed.

In [3]:

MultiSourceDataStore.display_config("scs2_config.yml")

MultiSourceDataStore.display_config("scs2_config.yml")

Configuration

User-defined ID	Data Store ID	Data Store Params	Data ID	Open Data Params	Grid-Mapping	Format
clms	clms	credentials: clms-credentials.json	tree-cover-density-2015|TCD_2015_020m_eu_03035_d05_E30N20.zarr	-	-	Zarr
clms	clms	credentials: clms-credentials.json	dominant-leaf-type-2015|DLT_2015_020m_eu_03035_d04_E30N20.zarr	-	-	Zarr
clms	clms	credentials: clms-credentials.json	forest-type-2015|FTY_2015_020m_eu_03035_d04_E30N20.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/annual_disturbances_1985_2023_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/disturbance_agent_1985_2023_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/greatest_disturbance_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/forest_mask_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/disturbance_agent_aggregated_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/latest_disturbance_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/number_disturbances_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/disturbance_probability_1985_2023_france.zarr	-	-	Zarr
senf	zenodo	root: 13333034	france/disturbance_severity_1985_2023_france.zarr	-	-	Zarr
gami	s3	root: dog.atlaseo-glm.eo-gridded-data; storage_options: {'anon': True, 'endpoint_url': 'https://s3.gfz-potsdam.de'}	collections/GAMI/GAMI_v2.1.zarr	-	-	Zarr
liu	zenodo	root: 8154445	planet_agb_30m_v0.1.tif	-	-	Zarr
liu	zenodo	root: 8154445	planet_canopy_cover_30m_v0.1.tif	-	-	Zarr
liu	zenodo	root: 8154445	planet_canopy_height_30m_v0.1.tif	-	-	Zarr

Now, we can initialize the MultiSourceDataStore by passing the path to the scs2_config.yml which currently is on the same level as this notebook.

By running the cell below, you would start seeing progress tables for each data that you requested in the scs2_config.yml.

NOTE: In the scs2_config.yml we are also using the custom_processing feature of this tool that allows us to run a function for processing each dataset separately. In this example, we have defined a module called modify_dataset that does some custom processing which takes a xarray.Dataset argument as input and returns a new xarray.Dataset object. To read more about this custom_processing function, you can see more here.

In [4]:

msds = MultiSourceDataStore("scs2_config.yml")

msds = MultiSourceDataStore("scs2_config.yml")

Preload Datasets from store 'zenodo_senf'

Dataset identifier	Status	Message	Exception
france.zip	COMPLETED	Already preloaded.	-

Preload Datasets from store 'clms'

Dataset identifier	Status	Message	Exception
tree-cover-density-2015|TCD_2015_020m_eu_03035_d05_E30N20	COMPLETED	Already preloaded.	-
dominant-leaf-type-2015|DLT_2015_020m_eu_03035_d04_E30N20	COMPLETED	Already preloaded.	-
forest-type-2015|FTY_2015_020m_eu_03035_d04_E30N20	COMPLETED	Already preloaded.	-

Cube Generation

Dataset identifier	Status	Message	Exception
clms	COMPLETED	Dataset 'clms' already generated.	-
senf	COMPLETED	Dataset 'senf' finished.	-
gami	COMPLETED	Dataset 'gami' already generated.	-
liu	COMPLETED	Dataset 'liu' already generated.	-

We can now open the data using the xcube datastore framework API as usual. Note that the multi-source data store requires a data store called storage, which is configured in our scs2_config.yml under the data_stores section.

In [5]:

ds = msds.stores.storage.open_data("gami.zarr")
ds

ds = msds.stores.storage.open_data("gami.zarr") ds

Out[5]:

<xarray.Dataset> Size: 3GB
Dimensions:      (members: 20, time: 2, latitude: 5199, longitude: 7000)
Coordinates:
  * latitude     (latitude) float64 42kB 44.8 44.8 44.8 44.8 ... 43.5 43.5 43.5
  * longitude    (longitude) float64 56kB -1.25 -1.25 -1.249 ... 0.4996 0.4999
  * members      (members) int64 160B 0 1 2 3 4 5 6 7 ... 13 14 15 16 17 18 19
    spatial_ref  int64 8B ...
  * time         (time) datetime64[ns] 16B 2010-01-01 2020-01-01
Data variables:
    forest_age   (members, time, latitude, longitude) int16 3GB dask.array<chunksize=(5, 1, 1000, 1000), meta=np.ndarray>
Attributes: (12/21)
    _FillValue:                 -9999
    contact:                    Simon Besnard (besnard@gfz.de) and Nuno Carva...
    created_by:                 Simon Besnard
    creation_date:              2025-01-31 16:52
    date_modified:              2025-09-11T14:51:07.671612
    frequency:                  2010 and 2020
    ...                         ...
    geospatial_lon_units:       degrees_east
    institute_id:               GFZ-Potsdam and MPI-BGC
    institution:                Helmholtz Center Potsdam GFZ German Research ...
    product_name:               Global Age Mapping Integration (GAMI) v2.1
    references:                 https://doi.org/10.5880/GFZ.1.4.2023.006 and ...
    version:                    v2.1

xarray.Dataset

• Dimensions:
  • members: 20
  • time: 2
  • latitude: 5199
  • longitude: 7000
• Coordinates: (5)
  • latitude
    (latitude)
    float64
    44.8 44.8 44.8 ... 43.5 43.5 43.5

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    array([44.799625, 44.799375, 44.799125, ..., 43.500625, 43.500375, 43.500125],
          shape=(5199,))

  • longitude
    (longitude)
    float64
    -1.25 -1.25 ... 0.4996 0.4999

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    array([-1.249875, -1.249625, -1.249375, ...,  0.499375,  0.499625,  0.499875],
          shape=(7000,))

  • members
    (members)
    int64
    0 1 2 3 4 5 6 ... 14 15 16 17 18 19

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19])

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    geographic_crs_name :

    WGS 84

    grid_mapping_name :

    latitude_longitude

    horizontal_datum_name :

    World Geodetic System 1984

    inverse_flattening :

    298.257223563

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    reference_ellipsoid_name :

    WGS 84

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    [1 values with dtype=int64]

  • time
    (time)
    datetime64[ns]
    2010-01-01 2020-01-01

    array(['2010-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • forest_age
    (members, time, latitude, longitude)
    int16
    dask.array<chunksize=(5, 1, 1000, 1000), meta=np.ndarray>

    grid_mapping :

    spatial_ref

    long_name :

    Forest age using a fusion machine learning and Landsat-based last time since disturbance

    units :

    years

    valid_max :

    300

    valid_min :

    1.0

    Array Chunk Bytes 2.71 GiB 9.54 MiB Shape (20, 2, 5199, 7000) (5, 1, 1000, 1000) Dask graph 336 chunks in 2 graph layers Data type int16 numpy.ndarray

• Indexes: (4)
  • latitude
    PandasIndex

    PandasIndex(Index([44.799625000000006, 44.799375000000005, 44.799125000000004,
            44.79887500000001,  44.79862500000001,  44.79837500000001,
           44.798125000000006, 44.797875000000005, 44.797625000000004,
            44.79737500000001,
           ...
           43.502374999999994,          43.502125,          43.501875,
                    43.501625, 43.501374999999996, 43.501124999999995,
            43.50087499999999,          43.500625,          43.500375,
                    43.500125],
          dtype='float64', name='latitude', length=5199))

  • longitude
    PandasIndex

    PandasIndex(Index([          -1.249875,           -1.249625, -1.2493750000000001,
                     -1.249125, -1.2488750000000002,           -1.248625,
                     -1.248375, -1.2481250000000002,           -1.247875,
                     -1.247625,
           ...
                      0.497625, 0.49787500000000007, 0.49812499999999993,
                      0.498375,  0.4986249999999999, 0.49887499999999996,
           0.49912500000000004,  0.4993749999999999,            0.499625,
                      0.499875],
          dtype='float64', name='longitude', length=7000))

  • members
    PandasIndex

    PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], dtype='int64', name='members'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2010-01-01', '2020-01-01'], dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (21)

  _FillValue :

  -9999

  contact :

  Simon Besnard (besnard@gfz.de) and Nuno Carvalhais (ncarvalhais@bgc-jena.mpg.de)

  created_by :

  Simon Besnard

  creation_date :

  2025-01-31 16:52

  date_modified :

  2025-09-11T14:51:07.671612

  frequency :

  2010 and 2020

  geospatial_bounds :

  POLYGON((-1.2524444450000096 43.49866666611111, -1.2524444450000096 44.802666667222226, 0.5022222227777894 44.802666667222226, 0.5022222227777894 43.49866666611111, -1.2524444450000096 43.49866666611111))

  geospatial_bounds_crs :

  CRS84

  geospatial_lat_max :

  44.802666667222226

  geospatial_lat_min :

  43.49866666611111

  geospatial_lat_resolution :

  0.00088889

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  0.5022222227777894

  geospatial_lon_min :

  -1.2524444450000096

  geospatial_lon_resolution :

  0.00088889

  geospatial_lon_units :

  degrees_east

  institute_id :

  GFZ-Potsdam and MPI-BGC

  institution :

  Helmholtz Center Potsdam GFZ German Research Centre for Geosciences and Max Planck Institute for Biogeochemistry

  product_name :

  Global Age Mapping Integration (GAMI) v2.1

  references :

  https://doi.org/10.5880/GFZ.1.4.2023.006 and https://doi.org/10.5194/essd-13-4881-2021

  version :

  v2.1

We can now select a variable for one timestep, downsample and plot it for a quick preview of the data

In [6]:

ds["forest_age"].isel(members=0, time=0)[::5, ::5].plot(vmin=0, vmax=120)

ds["forest_age"].isel(members=0, time=0)[::5, ::5].plot(vmin=0, vmax=120)

Out[6]:

<matplotlib.collections.QuadMesh at 0x7e1b4070f230>

In [7]:

ds["forest_age"].isel(members=0, time=1)[::5, ::5].plot(vmin=0, vmax=120)

ds["forest_age"].isel(members=0, time=1)[::5, ::5].plot(vmin=0, vmax=120)

Out[7]:

<matplotlib.collections.QuadMesh at 0x7e1b33d39450>

In [8]:

ds = msds.stores.storage.open_data("clms.zarr")
ds

ds = msds.stores.storage.open_data("clms.zarr") ds

Out[8]:

<xarray.Dataset> Size: 874MB
Dimensions:             (lat: 5199, lon: 7000)
Coordinates:
  * lat                 (lat) float64 42kB 44.8 44.8 44.8 ... 43.5 43.5 43.5
  * lon                 (lon) float64 56kB -1.25 -1.25 -1.249 ... 0.4996 0.4999
Data variables:
    dominant_leaf_type  (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    forest_type         (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    tree_cover_density  (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
Attributes:
    date_modified:              2025-09-11T15:03:36.475872
    geospatial_bounds:          POLYGON((-1.2167903401501168 43.3252872185624...
    geospatial_bounds_crs:      CRS84
    geospatial_lat_max:         44.972014235360945
    geospatial_lat_min:         43.32528721856248
    geospatial_lat_resolution:  0.00020337497095823664
    geospatial_lat_units:       degrees_north
    geospatial_lon_max:         0.47105574969819175
    geospatial_lon_min:         -1.2167903401501168
    geospatial_lon_resolution:  0.0002149804890727025
    geospatial_lon_units:       degrees_east

xarray.Dataset

• Dimensions:
  • lat: 5199
  • lon: 7000
• Coordinates: (2)
  • lat
    (lat)
    float64
    44.8 44.8 44.8 ... 43.5 43.5 43.5

    bounds :

    lat_bnds

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    array([44.799625, 44.799375, 44.799125, ..., 43.500625, 43.500375, 43.500125],
          shape=(5199,))

  • lon
    (lon)
    float64
    -1.25 -1.25 ... 0.4996 0.4999

    bounds :

    lon_bnds

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    array([-1.249875, -1.249625, -1.249375, ...,  0.499375,  0.499625,  0.499875],
          shape=(7000,))

• Data variables: (3)
  • dominant_leaf_type
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    DataType :

    Thematic

    RepresentationType :

    THEMATIC

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • forest_type
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    DataType :

    Thematic

    RepresentationType :

    THEMATIC

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • tree_cover_density
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    DataType :

    Thematic

    RepresentationType :

    THEMATIC

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([44.799625000000006, 44.799375000000005, 44.799125000000004,
            44.79887500000001,  44.79862500000001,  44.79837500000001,
           44.798125000000006, 44.797875000000005, 44.797625000000004,
            44.79737500000001,
           ...
           43.502374999999994,          43.502125,          43.501875,
                    43.501625, 43.501374999999996, 43.501124999999995,
            43.50087499999999,          43.500625,          43.500375,
                    43.500125],
          dtype='float64', name='lat', length=5199))

  • lon
    PandasIndex

    PandasIndex(Index([          -1.249875,           -1.249625, -1.2493750000000001,
                     -1.249125, -1.2488750000000002,           -1.248625,
                     -1.248375, -1.2481250000000002,           -1.247875,
                     -1.247625,
           ...
                      0.497625, 0.49787500000000007, 0.49812499999999993,
                      0.498375,  0.4986249999999999, 0.49887499999999996,
           0.49912500000000004,  0.4993749999999999,            0.499625,
                      0.499875],
          dtype='float64', name='lon', length=7000))

• Attributes: (11)

  date_modified :

  2025-09-11T15:03:36.475872

  geospatial_bounds :

  POLYGON((-1.2167903401501168 43.32528721856248, -1.2167903401501168 44.972014235360945, 0.47105574969819175 44.972014235360945, 0.47105574969819175 43.32528721856248, -1.2167903401501168 43.32528721856248))

  geospatial_bounds_crs :

  CRS84

  geospatial_lat_max :

  44.972014235360945

  geospatial_lat_min :

  43.32528721856248

  geospatial_lat_resolution :

  0.00020337497095823664

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  0.47105574969819175

  geospatial_lon_min :

  -1.2167903401501168

  geospatial_lon_resolution :

  0.0002149804890727025

  geospatial_lon_units :

  degrees_east

In [9]:

ds["dominant_leaf_type"][::5, ::5].plot(vmax=3)

ds["dominant_leaf_type"][::5, ::5].plot(vmax=3)

Out[9]:

<matplotlib.collections.QuadMesh at 0x7e1b34b6fed0>

In [10]:

ds["forest_type"][::5, ::5].plot(vmax=3)

ds["forest_type"][::5, ::5].plot(vmax=3)

Out[10]:

<matplotlib.collections.QuadMesh at 0x7e1b3d01a210>

In [11]:

ds["tree_cover_density"][::5, ::5].plot(vmax=100)

ds["tree_cover_density"][::5, ::5].plot(vmax=100)

Out[11]:

<matplotlib.collections.QuadMesh at 0x7e1b3d67c2d0>

In [12]:

ds = msds.stores.storage.open_data("liu.zarr")
ds

ds = msds.stores.storage.open_data("liu.zarr") ds

Out[12]:

<xarray.Dataset> Size: 874MB
Dimensions:        (lat: 5199, lon: 7000)
Coordinates:
  * lat            (lat) float64 42kB 44.8 44.8 44.8 44.8 ... 43.5 43.5 43.5
  * lon            (lon) float64 56kB -1.25 -1.25 -1.249 ... 0.4996 0.4999
Data variables:
    agb            (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    canopy_cover   (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    canopy_height  (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
Attributes:
    date_modified:              2025-09-11T14:52:36.048982
    geospatial_bounds:          POLYGON((-18.9594916134861 50.18279655877187,...
    geospatial_bounds_crs:      CRS84
    geospatial_lat_max:         52.03426577988996
    geospatial_lat_min:         50.18279655877187
    geospatial_lat_resolution:  0.00034892335219183224
    geospatial_lat_units:       degrees_north
    geospatial_lon_max:         -17.63688150306111
    geospatial_lon_min:         -18.9594916134861
    geospatial_lon_resolution:  0.0002391080323462802
    geospatial_lon_units:       degrees_east

xarray.Dataset

• Dimensions:
  • lat: 5199
  • lon: 7000
• Coordinates: (2)
  • lat
    (lat)
    float64
    44.8 44.8 44.8 ... 43.5 43.5 43.5

    bounds :

    lat_bnds

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    array([44.799625, 44.799375, 44.799125, ..., 43.500625, 43.500375, 43.500125],
          shape=(5199,))

  • lon
    (lon)
    float64
    -1.25 -1.25 ... 0.4996 0.4999

    bounds :

    lon_bnds

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    array([-1.249875, -1.249625, -1.249375, ...,  0.499375,  0.499625,  0.499875],
          shape=(7000,))

• Data variables: (3)
  • agb
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • canopy_cover
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • canopy_height
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([44.799625000000006, 44.799375000000005, 44.799125000000004,
            44.79887500000001,  44.79862500000001,  44.79837500000001,
           44.798125000000006, 44.797875000000005, 44.797625000000004,
            44.79737500000001,
           ...
           43.502374999999994,          43.502125,          43.501875,
                    43.501625, 43.501374999999996, 43.501124999999995,
            43.50087499999999,          43.500625,          43.500375,
                    43.500125],
          dtype='float64', name='lat', length=5199))

  • lon
    PandasIndex

    PandasIndex(Index([          -1.249875,           -1.249625, -1.2493750000000001,
                     -1.249125, -1.2488750000000002,           -1.248625,
                     -1.248375, -1.2481250000000002,           -1.247875,
                     -1.247625,
           ...
                      0.497625, 0.49787500000000007, 0.49812499999999993,
                      0.498375,  0.4986249999999999, 0.49887499999999996,
           0.49912500000000004,  0.4993749999999999,            0.499625,
                      0.499875],
          dtype='float64', name='lon', length=7000))

• Attributes: (11)

  date_modified :

  2025-09-11T14:52:36.048982

  geospatial_bounds :

  POLYGON((-18.9594916134861 50.18279655877187, -18.9594916134861 52.03426577988996, -17.63688150306111 52.03426577988996, -17.63688150306111 50.18279655877187, -18.9594916134861 50.18279655877187))

  geospatial_bounds_crs :

  CRS84

  geospatial_lat_max :

  52.03426577988996

  geospatial_lat_min :

  50.18279655877187

  geospatial_lat_resolution :

  0.00034892335219183224

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  -17.63688150306111

  geospatial_lon_min :

  -18.9594916134861

  geospatial_lon_resolution :

  0.0002391080323462802

  geospatial_lon_units :

  degrees_east

In [13]:

ds["agb"][::5, ::5].plot(vmax=20000)

ds["agb"][::5, ::5].plot(vmax=20000)

Out[13]:

<matplotlib.collections.QuadMesh at 0x7e1b3c645bd0>

In [14]:

ds["canopy_cover"][::5, ::5].plot()

ds["canopy_cover"][::5, ::5].plot()

Out[14]:

<matplotlib.collections.QuadMesh at 0x7e1b38ee8b90>

In [15]:

ds["canopy_height"][::5, ::5].plot()

ds["canopy_height"][::5, ::5].plot()

Out[15]:

<matplotlib.collections.QuadMesh at 0x7e1b430de350>

In [16]:

ds = msds.stores.storage.open_data("senf.zarr")
ds

ds = msds.stores.storage.open_data("senf.zarr") ds

Out[16]:

<xarray.Dataset> Size: 47GB
Dimensions:                       (time: 39, lat: 5199, lon: 7000)
Coordinates:
  * lat                           (lat) float64 42kB 44.8 44.8 ... 43.5 43.5
  * lon                           (lon) float64 56kB -1.25 -1.25 ... 0.4999
  * time                          (time) datetime64[ns] 312B 1985-01-01 ... 2...
Data variables:
    annual_disturbances           (time, lat, lon) float64 11GB dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
    disturbance_agent             (time, lat, lon) float64 11GB dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
    disturbance_agent_aggregated  (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    disturbance_probability       (time, lat, lon) float64 11GB dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
    disturbance_severity          (time, lat, lon) float64 11GB dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
    forest_mask                   (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    greatest_disturbance          (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    latest_disturbance            (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
    number_disturbances           (lat, lon) float64 291MB dask.array<chunksize=(1000, 1000), meta=np.ndarray>
Attributes:
    date_modified:              2025-09-11T15:19:16.815048
    geospatial_bounds:          POLYGON((-18.95915690343474 50.18273187952879...
    geospatial_bounds_crs:      CRS84
    geospatial_lat_max:         52.03419759001479
    geospatial_lat_min:         50.182731879528795
    geospatial_lat_resolution:  0.0003489227219688473
    geospatial_lat_units:       degrees_north
    geospatial_lon_max:         -17.6365344971688
    geospatial_lon_min:         -18.95915690343474
    geospatial_lon_resolution:  0.00023911036652535245
    geospatial_lon_units:       degrees_east

xarray.Dataset

• Dimensions:
  • time: 39
  • lat: 5199
  • lon: 7000
• Coordinates: (3)
  • lat
    (lat)
    float64
    44.8 44.8 44.8 ... 43.5 43.5 43.5

    bounds :

    lat_bnds

    long_name :

    latitude coordinate

    standard_name :

    latitude

    units :

    degrees_north

    array([44.799625, 44.799375, 44.799125, ..., 43.500625, 43.500375, 43.500125],
          shape=(5199,))

  • lon
    (lon)
    float64
    -1.25 -1.25 ... 0.4996 0.4999

    bounds :

    lon_bnds

    long_name :

    longitude coordinate

    standard_name :

    longitude

    units :

    degrees_east

    array([-1.249875, -1.249625, -1.249375, ...,  0.499375,  0.499625,  0.499875],
          shape=(7000,))

  • time
    (time)
    datetime64[ns]
    1985-01-01 ... 2023-01-01

    array(['1985-01-01T00:00:00.000000000', '1986-01-01T00:00:00.000000000',
           '1987-01-01T00:00:00.000000000', '1988-01-01T00:00:00.000000000',
           '1989-01-01T00:00:00.000000000', '1990-01-01T00:00:00.000000000',
           '1991-01-01T00:00:00.000000000', '1992-01-01T00:00:00.000000000',
           '1993-01-01T00:00:00.000000000', '1994-01-01T00:00:00.000000000',
           '1995-01-01T00:00:00.000000000', '1996-01-01T00:00:00.000000000',
           '1997-01-01T00:00:00.000000000', '1998-01-01T00:00:00.000000000',
           '1999-01-01T00:00:00.000000000', '2000-01-01T00:00:00.000000000',
           '2001-01-01T00:00:00.000000000', '2002-01-01T00:00:00.000000000',
           '2003-01-01T00:00:00.000000000', '2004-01-01T00:00:00.000000000',
           '2005-01-01T00:00:00.000000000', '2006-01-01T00:00:00.000000000',
           '2007-01-01T00:00:00.000000000', '2008-01-01T00:00:00.000000000',
           '2009-01-01T00:00:00.000000000', '2010-01-01T00:00:00.000000000',
           '2011-01-01T00:00:00.000000000', '2012-01-01T00:00:00.000000000',
           '2013-01-01T00:00:00.000000000', '2014-01-01T00:00:00.000000000',
           '2015-01-01T00:00:00.000000000', '2016-01-01T00:00:00.000000000',
           '2017-01-01T00:00:00.000000000', '2018-01-01T00:00:00.000000000',
           '2019-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000',
           '2021-01-01T00:00:00.000000000', '2022-01-01T00:00:00.000000000',
           '2023-01-01T00:00:00.000000000'], dtype='datetime64[ns]')

• Data variables: (9)
  • annual_disturbances
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>

    source :

    file:///home/konstantin/bc_kon/01_coding/01_github/xcube-multistore/examples/example2/zenodo_cache/13333034/downloads/france/annual_disturbances_1985_2023_france.tif

    Array Chunk Bytes 10.57 GiB 7.63 MiB Shape (39, 5199, 7000) (1, 1000, 1000) Dask graph 1638 chunks in 2 graph layers Data type float64 numpy.ndarray

  • disturbance_agent
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>

    source :

    file:///home/konstantin/bc_kon/01_coding/01_github/xcube-multistore/examples/example2/zenodo_cache/13333034/downloads/france/disturbance_agent_1985_2023_france.tif

    Array Chunk Bytes 10.57 GiB 7.63 MiB Shape (39, 5199, 7000) (1, 1000, 1000) Dask graph 1638 chunks in 2 graph layers Data type float64 numpy.ndarray

  • disturbance_agent_aggregated
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • disturbance_probability
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>

    source :

    file:///home/konstantin/bc_kon/01_coding/01_github/xcube-multistore/examples/example2/zenodo_cache/13333034/downloads/france/disturbance_probability_1985_2023_france.tif

    Array Chunk Bytes 10.57 GiB 7.63 MiB Shape (39, 5199, 7000) (1, 1000, 1000) Dask graph 1638 chunks in 2 graph layers Data type float64 numpy.ndarray

  • disturbance_severity
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>

    source :

    file:///home/konstantin/bc_kon/01_coding/01_github/xcube-multistore/examples/example2/zenodo_cache/13333034/downloads/france/disturbance_severity_1985_2023_france.tif

    Array Chunk Bytes 10.57 GiB 7.63 MiB Shape (39, 5199, 7000) (1, 1000, 1000) Dask graph 1638 chunks in 2 graph layers Data type float64 numpy.ndarray

  • forest_mask
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • greatest_disturbance
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • latest_disturbance
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

  • number_disturbances
    (lat, lon)
    float64
    dask.array<chunksize=(1000, 1000), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    Array Chunk Bytes 277.66 MiB 7.63 MiB Shape (5199, 7000) (1000, 1000) Dask graph 42 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([44.799625000000006, 44.799375000000005, 44.799125000000004,
            44.79887500000001,  44.79862500000001,  44.79837500000001,
           44.798125000000006, 44.797875000000005, 44.797625000000004,
            44.79737500000001,
           ...
           43.502374999999994,          43.502125,          43.501875,
                    43.501625, 43.501374999999996, 43.501124999999995,
            43.50087499999999,          43.500625,          43.500375,
                    43.500125],
          dtype='float64', name='lat', length=5199))

  • lon
    PandasIndex

    PandasIndex(Index([          -1.249875,           -1.249625, -1.2493750000000001,
                     -1.249125, -1.2488750000000002,           -1.248625,
                     -1.248375, -1.2481250000000002,           -1.247875,
                     -1.247625,
           ...
                      0.497625, 0.49787500000000007, 0.49812499999999993,
                      0.498375,  0.4986249999999999, 0.49887499999999996,
           0.49912500000000004,  0.4993749999999999,            0.499625,
                      0.499875],
          dtype='float64', name='lon', length=7000))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1985-01-01', '1986-01-01', '1987-01-01', '1988-01-01',
                   '1989-01-01', '1990-01-01', '1991-01-01', '1992-01-01',
                   '1993-01-01', '1994-01-01', '1995-01-01', '1996-01-01',
                   '1997-01-01', '1998-01-01', '1999-01-01', '2000-01-01',
                   '2001-01-01', '2002-01-01', '2003-01-01', '2004-01-01',
                   '2005-01-01', '2006-01-01', '2007-01-01', '2008-01-01',
                   '2009-01-01', '2010-01-01', '2011-01-01', '2012-01-01',
                   '2013-01-01', '2014-01-01', '2015-01-01', '2016-01-01',
                   '2017-01-01', '2018-01-01', '2019-01-01', '2020-01-01',
                   '2021-01-01', '2022-01-01', '2023-01-01'],
                  dtype='datetime64[ns]', name='time', freq=None))

• Attributes: (11)

  date_modified :

  2025-09-11T15:19:16.815048

  geospatial_bounds :

  POLYGON((-18.95915690343474 50.182731879528795, -18.95915690343474 52.03419759001479, -17.6365344971688 52.03419759001479, -17.6365344971688 50.182731879528795, -18.95915690343474 50.182731879528795))

  geospatial_bounds_crs :

  CRS84

  geospatial_lat_max :

  52.03419759001479

  geospatial_lat_min :

  50.182731879528795

  geospatial_lat_resolution :

  0.0003489227219688473

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  -17.6365344971688

  geospatial_lon_min :

  -18.95915690343474

  geospatial_lon_resolution :

  0.00023911036652535245

  geospatial_lon_units :

  degrees_east

In [17]:

ds["annual_disturbances"].isel(time=10)[::5, ::5].plot()

ds["annual_disturbances"].isel(time=10)[::5, ::5].plot()

Out[17]:

<matplotlib.collections.QuadMesh at 0x7e1b3029ca50>

In [18]:

ds["disturbance_probability"].isel(time=10)[::5, ::5].plot()

ds["disturbance_probability"].isel(time=10)[::5, ::5].plot()

Out[18]:

<matplotlib.collections.QuadMesh at 0x7e1b40412210>