Load Copernicus Land Monitoring Service (CLMS) data

from odc import stac as odc_stac
from pystac_client import Client
import xarray as xr
from pyproj import Transformer
import math
import xarray as xr
import numpy as np
from rasterio.windows import from_bounds

Connect to eodc stac and open the collection and the item for the Copernicus Land Monitoring data

eodc_catalog=Client.open("https://stac.eodc.eu/api/v1")
col = eodc_catalog.get_collection("clms")
col

• type "Collection"
• id "clms"
• stac_version "1.1.0"
• description "The Copernicus Land Monitoring Service (CLMS) is a component of the European Union’s Copernicus Earth observation programme that provides systematic information on land cover, land use, and biophysical parameters across Europe and globally. It delivers consistent, long-term datasets derived primarily from satellite imagery to support environmental monitoring, spatial planning, agriculture, forestry, and climate applications. The service produces several product lines, including pan-European, local, and global components, enabling both high-resolution mapping and large-scale assessments."
• links[] 5 items

  • 0

    • rel "items"
    • href "https://stac.eodc.eu/api/v1/collections/clms/items"
    • type "application/geo+json"

  • 1

    • rel "parent"
    • href "https://stac.eodc.eu/api/v1/"
    • type "application/json"

  • 2

    • rel "root"
    • href "https://stac.eodc.eu/api/v1"
    • type "application/json"
    • title "EODC Data Catalogue"

  • 3

    • rel "self"
    • href "https://stac.eodc.eu/api/v1/collections/clms"
    • type "application/json"

  • 4

    • rel "http://www.opengis.net/def/rel/ogc/1.0/queryables"
    • href "https://stac.eodc.eu/api/v1/collections/clms/queryables"
    • type "application/schema+json"
    • title "Queryables"
• stac_extensions[] 1 items

  • 0 "https://stac-extensions.github.io/zarr/v1.1.0/schema.json"
• title "Copernicus Land Monitoring Service"
• extent

  • spatial

    • bbox[] 1 items

      • 0[] 4 items

        • 0 8.45

        • 1 45.6

        • 2 18.07

        • 3 49.83
  • temporal

    • interval[] 1 items

      • 0[] 2 items

        • 0 "2019-01-01T00:00:00Z"

        • 1 "2022-12-31T00:00:00Z"
• license "CC-BY-4.0"
• assets

  • thumbnail

    • href "https://raw.githubusercontent.com/koenifra/CLMS_zarr/main/CLMS_thumbnail.png"
    • type "image/png"
    • title "CLMS thumbnail"
    • roles[] 1 items

      • 0 "thumbnail"

search = eodc_catalog.search(collections=["clms"])
for item in search.items():
    print(item.id)

clms-vpp
clms-st

item = next(
    eodc_catalog.search(collections=["clms"], ids=["clms-vpp"]).get_items()
)

item

/home/katharina/.pyenv/versions/3.11.5/lib/python3.11/site-packages/pystac_client/item_search.py:925: FutureWarning: get_items() is deprecated, use items() instead
  warnings.warn(

• type "Feature"
• stac_version "1.1.0"
• stac_extensions[] 0 items
• id "clms-vpp"
• geometry

  • type "Polygon"
  • coordinates[] 1 items

    • 0[] 5 items

      • 0[] 2 items

        • 0 8.45

        • 1 45.6

      • 1[] 2 items

        • 0 18.07

        • 1 45.6

      • 2[] 2 items

        • 0 18.07

        • 1 49.83

      • 3[] 2 items

        • 0 8.45

        • 1 49.83

      • 4[] 2 items

        • 0 8.45

        • 1 45.6
• bbox[] 4 items

  • 0 8.45

  • 1 45.6

  • 2 18.07

  • 3 49.83
• properties

  • title "Vegetation Phenology and Productivity"
  • description "HR-VPP product suite provides information on vegetation status and dynamics."
  • end_datetime "2022-12-31T00:00:00Z"
  • start_datetime "2019-01-01T00:00:00Z"
  • datetime None
• links[] 4 items

  • 0

    • rel "collection"
    • href "https://stac.eodc.eu/api/v1/collections/clms"
    • type "application/json"

  • 1

    • rel "parent"
    • href "https://stac.eodc.eu/api/v1/collections/clms"
    • type "application/json"

  • 2

    • rel "root"
    • href "https://stac.eodc.eu/api/v1"
    • type "application/json"
    • title "EODC Data Catalogue"

  • 3

    • rel "self"
    • href "https://stac.eodc.eu/api/v1/collections/clms/items/clms-vpp"
    • type "application/geo+json"
• assets

  • data

    • href "https://data.eodc.eu/collections/CLMS/CLMS.zarr"
    • type "application/vnd.zarr; version=3"
    • title "Vegetation Phenology and Productivity"
    • description "HR-VPP product suite provides information on vegetation status and dynamics."
    • bands[] 13 items

      • 0

        • name "AMPL"
        • description "Season Amplitude"

      • 1

        • name "EOSD"
        • description "End-of-season Date"

      • 2

        • name "EOSV"
        • description "End-of-season Value"

      • 3

        • name "LENGTH"
        • description "Season Length"

      • 4

        • name "LSLOPE"
        • description "Slope of the Green-up Period"

      • 5

        • name "MAXV"
        • description "Season Maximum Value"

      • 6

        • name "MINV"
        • description "Season Minimum Value"

      • 7

        • name "QFLAG"
        • description "Quality Flag"

      • 8

        • name "RSLOPE"
        • description "Slope of the Green-down Period"

      • 9

        • name "SOSD"
        • description "Start-of-season Date"

      • 10

        • name "SOSV"
        • description "Start-of-season Value"

      • 11

        • name "SPROD"
        • description "Seasonal Productivity"

      • 12

        • name "TPROD"
        • description "Total Productivity"
    • xarray:open_zarr_kwargs

      • group "VPP"
      • consolidated True
      • zarr_version 3
    • roles[] 1 items

      • 0 "data"
• collection "clms"

asset = item.assets["data"]
asset

• href "https://data.eodc.eu/collections/CLMS/CLMS.zarr"
• type "application/vnd.zarr; version=3"
• title "Vegetation Phenology and Productivity"
• description "HR-VPP product suite provides information on vegetation status and dynamics."
• bands[] 13 items

  • 0

    • name "AMPL"
    • description "Season Amplitude"

  • 1

    • name "EOSD"
    • description "End-of-season Date"

  • 2

    • name "EOSV"
    • description "End-of-season Value"

  • 3

    • name "LENGTH"
    • description "Season Length"

  • 4

    • name "LSLOPE"
    • description "Slope of the Green-up Period"

  • 5

    • name "MAXV"
    • description "Season Maximum Value"

  • 6

    • name "MINV"
    • description "Season Minimum Value"

  • 7

    • name "QFLAG"
    • description "Quality Flag"

  • 8

    • name "RSLOPE"
    • description "Slope of the Green-down Period"

  • 9

    • name "SOSD"
    • description "Start-of-season Date"

  • 10

    • name "SOSV"
    • description "Start-of-season Value"

  • 11

    • name "SPROD"
    • description "Seasonal Productivity"

  • 12

    • name "TPROD"
    • description "Total Productivity"
• xarray:open_zarr_kwargs

  • group "VPP"
  • consolidated True
  • zarr_version 3
• roles[] 1 items

  • 0 "data"

print(list(item.assets.keys()))
asset = item.assets["data"]
print(asset.href)

['data']
https://data.eodc.eu/collections/CLMS/CLMS.zarr

store_url = asset.href
group_path = asset.extra_fields.get("zarr:group")  # may be None

print("href:", store_url)
print("group:", group_path)

href: https://data.eodc.eu/collections/CLMS/CLMS.zarr
group: None

Open the zarr file

coords = xr.open_zarr(
    store_url,
    **asset.extra_fields["xarray:open_zarr_kwargs"]
)

/tmp/ipykernel_1935098/348867924.py:1: FutureWarning: zarr_version is deprecated, use zarr_format
  coords = xr.open_zarr(

maxv_ds = xr.open_zarr(
    store_url,
    group="VPP/MAXV",
    zarr_version=3,
    chunks="auto",
)

/tmp/ipykernel_1935098/2837227757.py:1: FutureWarning: zarr_version is deprecated, use zarr_format
  maxv_ds = xr.open_zarr(

maxv_zarr = maxv_ds["SEASON1"].assign_coords(
    time=("time", coords["time"].values),
    x=("x", coords["x"].values),
    y=("y", coords["y"].values),
)

Select an area: Here Rosalia

The data is in EPSG3035, so we need to transform the corrdinates.

minlon, minlat, maxlon, maxlat = 16.25, 46.68, 16.31, 47.72
year = 2022

transformer = Transformer.from_crs(
    "EPSG:4326",
    "EPSG:3035",
    always_xy=True
)

def to_3035(minlon, minlat, maxlon, maxlat):
    x0, y0 = transformer.transform(minlon, minlat)
    x1, y1 = transformer.transform(maxlon, minlat)
    x2, y2 = transformer.transform(maxlon, maxlat)
    x3, y3 = transformer.transform(minlon, maxlat)

    minx_3035 = int(math.floor(min(x0, x1, x2, x3)))
    maxx_3035 = int(math.ceil (max(x0, x1, x2, x3)))
    miny_3035 = int(math.floor(min(y0, y1, y2, y3)))
    maxy_3035 = int(math.ceil (max(y0, y1, y2, y3)))

    return minx_3035, miny_3035, maxx_3035, maxy_3035

minx, miny, maxx, maxy = to_3035(minlon, minlat, maxlon, maxlat)

print(minx, miny, maxx, maxy)

4789576 2638650 4803578 2754238

Get the data in the selected area

subset = maxv_zarr.sel(
    time=year,
    x=slice(minx, maxx),
    y=slice(miny, maxy),
)

print(subset.values)

[[0.6665 0.6759 1.0142 ... 1.61   1.554  1.7099]
 [0.5018 1.0029 1.3698 ... 2.0067 1.8184 1.4919]
 [0.7442 1.3635 1.6927 ... 2.1152 2.3205 1.9617]
 ...
 [1.4486 2.0843 2.2482 ... 1.2622 1.2122 1.3095]
 [1.3237 1.3887 1.3536 ... 1.2472 1.0891 0.9634]
 [1.1522 1.0731 0.9404 ... 0.9311 0.806  0.6154]]