Access data from nonsearchable stac catalog

Load a GeoTiff file from a non-seachable STAC catalog¶

A DeepESDL example notebook¶

This notebook demonstrates how to load a GeoTIFF file from a non-searchable STAC catalog via xcube-stac store. The data in the example is fetchedthe EcoDataCube.eu STAC catalog.

A non-searchable catalog does not implement the STAC API - Item Search conformance class. When searching in such type of catalog, the catalog needs to be crawled through and the items properties needs to be matched to the search parameters. This process can be slow, especially for large catalogs.

For more example notebooks using xcube-stac, please head over to GitHub https://github.com/xcube-dev/xcube-stac/tree/main/examples/notebooks. There you can find examples like leading GeoTiffs and netCDFs from searchable STAC-Catalogs.

Please, also refer to the DeepESDL documentation and visit the platform's website for further information!

Brockmann Consult, 2025

---

This notebook runs with the python environment users-deepesdl-xcube-1.11.1, please checkout the documentation for help on changing the environment.

In [1]:

%%time
from xcube.core.store import new_data_store, get_data_store_params_schema
import itertools

%%time from xcube.core.store import new_data_store, get_data_store_params_schema import itertools

CPU times: user 1.5 s, sys: 519 ms, total: 2.02 s
Wall time: 3.94 s

First, we get the store parameters needed to initialize a STAC data store.

In [2]:

get_data_store_params_schema("stac")

get_data_store_params_schema("stac")

Out[2]:

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

We determine the url of the EcoDataCube.eu STAC catalog and initiate a STAC data store where the xcube-stac plugin is recognized by setting the first argument to "stac" in the new_data_store function.

In [3]:

%%time
url = "https://s3.eu-central-1.wasabisys.com/stac/odse/catalog.json"
store = new_data_store("stac", url=url)

%%time url = "https://s3.eu-central-1.wasabisys.com/stac/odse/catalog.json" store = new_data_store("stac", url=url)

CPU times: user 29.1 ms, sys: 1.59 ms, total: 30.7 ms
Wall time: 262 ms

/home/conda/users/60ba36d6-1759325445-137-deepesdl-xcube-1.11.1/lib/python3.13/site-packages/pystac_client/client.py:191: NoConformsTo: Server does not advertise any conformance classes.
  warnings.warn(NoConformsTo())

The data IDs point to a STAC item's JSON and are specified by the segment of the URL that follows the catalog's URL. The data IDs can be streamed using the following code where we show the first 10 data IDs as an example.

In [4]:

%%time
data_ids = store.get_data_ids()
list(itertools.islice(data_ids, 10))

%%time data_ids = store.get_data_ids() list(itertools.islice(data_ids, 10))

CPU times: user 1.63 s, sys: 97.9 ms, total: 1.72 s
Wall time: 25.6 s

Out[4]:

['accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000101_20000131/accum.precipitation_chelsa.montlhy_20000101_20000131.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000201_20000228/accum.precipitation_chelsa.montlhy_20000201_20000228.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000301_20000331/accum.precipitation_chelsa.montlhy_20000301_20000331.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000401_20000430/accum.precipitation_chelsa.montlhy_20000401_20000430.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000501_20000531/accum.precipitation_chelsa.montlhy_20000501_20000531.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000601_20000630/accum.precipitation_chelsa.montlhy_20000601_20000630.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000701_20000731/accum.precipitation_chelsa.montlhy_20000701_20000731.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000801_20000831/accum.precipitation_chelsa.montlhy_20000801_20000831.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20000901_20000930/accum.precipitation_chelsa.montlhy_20000901_20000930.json',
 'accum.precipitation_chelsa.montlhy/accum.precipitation_chelsa.montlhy_20001001_20001031/accum.precipitation_chelsa.montlhy_20001001_20001031.json']

In the next step, we can search for items using search parameters. The following code shows which search parameters are available.

In [5]:

%%time
search_params = store.get_search_params_schema()
search_params

%%time search_params = store.get_search_params_schema() search_params

CPU times: user 23 μs, sys: 6 μs, total: 29 μs
Wall time: 31.7 μs

Out[5]:

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

Now, let's search for Topographic Wetness Index generated by the Ensemble Digital Terrain Model by Whitbox Workflow for the European region during the first quarter of 2010.

In [6]:

%%time
descriptors = list(
    store.search_data(
        collections=["twi_edtm"],
        bbox=[-10, 40, 40, 70],
        time_range=["2010-01-01", "2010-04-01"],
    )
)
[d.to_dict() for d in descriptors]

%%time descriptors = list( store.search_data( collections=["twi_edtm"], bbox=[-10, 40, 40, 70], time_range=["2010-01-01", "2010-04-01"], ) ) [d.to_dict() for d in descriptors]

CPU times: user 32.5 ms, sys: 0 ns, total: 32.5 ms
Wall time: 183 ms

Out[6]:

[{'data_id': 'twi_edtm/twi_edtm_20060101_20151231/twi_edtm_20060101_20151231.json',
  'data_type': 'dataset',
  'bbox': [-56.51881139294227,
   24.275788389340878,
   72.93153043046625,
   72.64259665547773],
  'time_range': ('2006-01-01', '2015-12-31')}]

In the next step, we can open the data for each data ID. The following code shows which parameters are available for opening the data.

In [7]:

%%time
open_params = store.get_open_data_params_schema()
open_params

%%time open_params = store.get_open_data_params_schema() open_params

CPU times: user 39 μs, sys: 0 ns, total: 39 μs
Wall time: 42.2 μs

Out[7]:

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

Now, we open the data for a given data ID, where we select all available assets, which can be opened by the data store. Note that we assign data_type to get an xarray.Dataset as return value.

In [8]:

%%time
ds = store.open_data(descriptors[0].data_id, data_type="dataset")
ds

%%time ds = store.open_data(descriptors[0].data_id, data_type="dataset") ds

CPU times: user 259 ms, sys: 71.6 ms, total: 331 ms
Wall time: 975 ms

Out[8]:

<xarray.Dataset> Size: 66GB
Dimensions:      (x: 216700, y: 153400)
Coordinates:
  * x            (x) float64 2MB 9e+05 9e+05 9.001e+05 ... 7.401e+06 7.401e+06
  * y            (y) float64 1MB 5.501e+06 5.501e+06 ... 8.99e+05 8.99e+05
    spatial_ref  int64 8B 0
Data variables:
    band_1       (y, x) int16 66GB dask.array<chunksize=(512, 512), meta=np.ndarray>
Attributes:
    stac_catalog_url:    https://s3.eu-central-1.wasabisys.com/stac/odse/cata...
    stac_item_id:        twi_edtm_20060101_20151231
    xcube_stac_version:  1.1.1

xarray.Dataset

• Dimensions:
  • x: 216700
  • y: 153400
• Coordinates: (3)
  • x
    (x)
    float64
    9e+05 9e+05 ... 7.401e+06 7.401e+06

    array([ 900015.,  900045.,  900075., ..., 7400925., 7400955., 7400985.],
          shape=(216700,))

  • y
    (y)
    float64
    5.501e+06 5.501e+06 ... 8.99e+05

    array([5500985., 5500955., 5500925., ...,  899075.,  899045.,  899015.],
          shape=(153400,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCRS["ETRS89-extended / LAEA Europe",BASEGEOGCRS["ETRS89",DATUM["European Terrestrial Reference System 1989",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4258]],CONVERSION["unnamed",METHOD["Lambert Azimuthal Equal Area",ID["EPSG",9820]],PARAMETER["Latitude of natural origin",52,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",10,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",4321000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",3210000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["northing",north,ORDER[1],LENGTHUNIT["metre",1]],AXIS["easting",east,ORDER[2],LENGTHUNIT["metre",1]],ID["EPSG",3035]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    ETRS89

    horizontal_datum_name :

    European Terrestrial Reference System 1989

    projected_crs_name :

    ETRS89-extended / LAEA Europe

    grid_mapping_name :

    lambert_azimuthal_equal_area

    latitude_of_projection_origin :

    52.0

    longitude_of_projection_origin :

    10.0

    false_easting :

    4321000.0

    false_northing :

    3210000.0

    array(0)

• Data variables: (1)
  • band_1
    (y, x)
    int16
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    32767

    scale_factor :

    1.0

    add_offset :

    0.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 61.92 GiB 512.00 kiB Shape (153400, 216700) (512, 512) Dask graph 127200 chunks in 3 graph layers Data type int16 numpy.ndarray

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([ 900015.0,  900045.0,  900075.0,  900105.0,  900135.0,  900165.0,
            900195.0,  900225.0,  900255.0,  900285.0,
           ...
           7400715.0, 7400745.0, 7400775.0, 7400805.0, 7400835.0, 7400865.0,
           7400895.0, 7400925.0, 7400955.0, 7400985.0],
          dtype='float64', name='x', length=216700))

  • y
    PandasIndex

    PandasIndex(Index([5500985.0, 5500955.0, 5500925.0, 5500895.0, 5500865.0, 5500835.0,
           5500805.0, 5500775.0, 5500745.0, 5500715.0,
           ...
            899285.0,  899255.0,  899225.0,  899195.0,  899165.0,  899135.0,
            899105.0,  899075.0,  899045.0,  899015.0],
          dtype='float64', name='y', length=153400))

• Attributes: (3)

  stac_catalog_url :

  https://s3.eu-central-1.wasabisys.com/stac/odse/catalog.json

  stac_item_id :

  twi_edtm_20060101_20151231

  xcube_stac_version :

  1.1.1

We plot the loaded data as an example below.

In [9]:

%%time
ds.band_1[100000:120000:20, 100000:120000:20].plot(vmin=0, vmax=800)

%%time ds.band_1[100000:120000:20, 100000:120000:20].plot(vmin=0, vmax=800)

CPU times: user 4.01 s, sys: 1.14 s, total: 5.16 s
Wall time: 16 s

Out[9]:

<matplotlib.collections.QuadMesh at 0x7f1ae9960830>

We can also open a GeoTIFF as a xcube's multi-resolution dataset, where we can select the level of resolution, shown below.

In [10]:

%%time
mlds = store.open_data(descriptors[0].data_id, data_type="mldataset")
mlds.num_levels

%%time mlds = store.open_data(descriptors[0].data_id, data_type="mldataset") mlds.num_levels

CPU times: user 37 ms, sys: 3.96 ms, total: 40.9 ms
Wall time: 277 ms

Out[10]:

8

In [11]:

ds = mlds.get_dataset(5)
ds

ds = mlds.get_dataset(5) ds

Out[11]:

<xarray.Dataset> Size: 65MB
Dimensions:      (x: 6771, y: 4793)
Coordinates:
  * x            (x) float64 54kB 9.005e+05 9.014e+05 ... 7.4e+06 7.401e+06
  * y            (y) float64 38kB 5.501e+06 5.5e+06 ... 9.004e+05 8.995e+05
    spatial_ref  int64 8B 0
Data variables:
    band_1       (y, x) int16 65MB dask.array<chunksize=(512, 512), meta=np.ndarray>
Attributes:
    source:   https://s3.ecodatacube.eu/arco/twi_edtm_m_30m_s_20000101_202212...

xarray.Dataset

• Dimensions:
  • x: 6771
  • y: 4793
• Coordinates: (3)
  • x
    (x)
    float64
    9.005e+05 9.014e+05 ... 7.401e+06

    array([ 900480.062029,  901440.186088,  902400.310146, ..., 7398599.689854,
           7399559.813912, 7400519.937971], shape=(6771,))

  • y
    (y)
    float64
    5.501e+06 5.5e+06 ... 8.995e+05

    array([5500519.92489 , 5499559.774671, 5498599.624452, ...,  901400.375548,
            900440.225329,  899480.07511 ], shape=(4793,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314140356

    inverse_flattening :

    298.257222101

    reference_ellipsoid_name :

    GRS 1980

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    ETRS89

    horizontal_datum_name :

    European Terrestrial Reference System 1989

    projected_crs_name :

    ETRS89-extended / LAEA Europe

    grid_mapping_name :

    lambert_azimuthal_equal_area

    latitude_of_projection_origin :

    52.0

    longitude_of_projection_origin :

    10.0

    false_easting :

    4321000.0

    false_northing :

    3210000.0

    spatial_ref :

    PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]

    GeoTransform :

    900000.0 960.1240584847143 0.0 5501000.0 0.0 -960.1502190694762

    array(0)

• Data variables: (1)
  • band_1
    (y, x)
    int16
    dask.array<chunksize=(512, 512), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    32767

    scale_factor :

    1.0

    add_offset :

    0.0

    grid_mapping :

    spatial_ref

    Array Chunk Bytes 61.90 MiB 512.00 kiB Shape (4793, 6771) (512, 512) Dask graph 140 chunks in 3 graph layers Data type int16 numpy.ndarray

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([900480.0620292424, 901440.1860877271, 902400.3101462119,
           903360.4342046966, 904320.5582631812, 905280.6823216659,
           906240.8063801507, 907200.9304386354, 908161.0544971201,
           909121.1785556048,
           ...
           7391878.821444396, 7392838.945502881, 7393799.069561365,
            7394759.19361985, 7395719.317678334, 7396679.441736819,
           7397639.565795304, 7398599.689853788, 7399559.813912273,
           7400519.937970758],
          dtype='float64', name='x', length=6771))

  • y
    PandasIndex

    PandasIndex(Index([5500519.924890465, 5499559.774671395, 5498599.624452326,
           5497639.474233257, 5496679.324014187, 5495719.173795118,
           5494759.023576048, 5493798.873356978,  5492838.72313791,
            5491878.57291884,
           ...
           908121.4270811602, 907161.2768620914, 906201.1266430216,
           905240.9764239518, 904280.8262048829, 903320.6759858131,
           902360.5257667433, 901400.3755476745, 900440.2253286047,
           899480.0751095349],
          dtype='float64', name='y', length=4793))

• Attributes: (1)

  source :

  https://s3.ecodatacube.eu/arco/twi_edtm_m_30m_s_20000101_20221231_eu_epsg.3035_v20241230.tif

In [12]:

%%time
ds.band_1[3125:3750, 3125:3750].plot(vmax=800)

%%time ds.band_1[3125:3750, 3125:3750].plot(vmax=800)

CPU times: user 85 ms, sys: 10.3 ms, total: 95.3 ms
Wall time: 377 ms

Out[12]:

<matplotlib.collections.QuadMesh at 0x7f1aed716d50>

In [ ]: