STAC connect to ClimateDT weather stations

from odc import stac as odc_stac
from pystac_client import Client
import xarray as xr
import xarray as xr
import numpy as np

Connect to the eodc STAC catalogue with pystac-client

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

collection = "climatedt-austria"

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

IFS-NEMO-ScenarioMIP-snowDepth
IFS-NEMO-ScenarioMIP
ICON-ScenarioMIP
sd-IFS-FESOM-story-nudging
IFS-FESOM-story-nudging-Tplus2
IFS-FESOM-story-nudging-hist
IFS-FESOM-story-nudging-cont
2t-IFS-FESOM-story-nudging
ICON-CMIP6
IFS-NEMO-HighResMIP

IFS-NEMO-CMIP6

item = next(
    eodc_catalog.search(collections=[collection], ids=["IFS-FESOM-story-nudging-Tplus2"]).items()
)
item
<Item id=IFS-FESOM-story-nudging-Tplus2> 
asset = item.assets["data"]
store_url = asset.href
print("href:", store_url)

href: https://objects.eodc.eu/68e13833a1624f43ba2cac01376a18af:destine-climate-dt/Austria/climate-dt-austria.zarr

Open data with xarray

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

print(coords)

<xarray.Dataset> Size: 508MB
Dimensions:    (t: 61344, pointid: 517)
Coordinates:
  * t          (t) datetime64[ns] 491kB 2017-01-01 ... 2023-12-31T23:00:00
  * pointid    (pointid) int64 4kB 1 0 4 9 6 5 7 ... 510 511 512 513 514 515 516
    latitude   (pointid) float64 4kB dask.array<chunksize=(517,), meta=np.ndarray>
    longitude  (pointid) float64 4kB dask.array<chunksize=(517,), meta=np.ndarray>
Data variables:
    2t         (t, pointid) float64 254MB dask.array<chunksize=(61344, 100), meta=np.ndarray>
    tp         (t, pointid) float64 254MB dask.array<chunksize=(61344, 100), meta=np.ndarray>

coords["2t"].isel(pointid=99).plot()

[<matplotlib.lines.Line2D at 0x76f157d79610>]
../_images/8f6f3d08c153459bc21fe945ee497ecee527ae1019675d1c1e6cfdb582f732bc.png

Create subset with bbox and year

# Bounding box
lon_min, lat_min, lon_max, lat_max = 15.839539,47.565407,16.877747,48.039529
# Compute lat/lon (small, so this is fast)
lats = coords.latitude.compute()
lons = coords.longitude.compute()

# Boolean mask for points within the bbox
mask = (
    (lats >= lat_min) & (lats <= lat_max) &
    (lons >= lon_min) & (lons <= lon_max)
)

# Subset
ds_subset = coords.isel(pointid=mask.values)

ds_2023 = ds_subset.sel(t=ds_subset.t.dt.year == 2023).compute()
print(ds_2023)

<xarray.Dataset> Size: 3MB
Dimensions:    (t: 8760, pointid: 19)
Coordinates:
  * t          (t) datetime64[ns] 70kB 2023-01-01 ... 2023-12-31T23:00:00
  * pointid    (pointid) int64 152B 310 311 312 328 330 ... 369 371 370 372 373
    latitude   (pointid) float64 152B 47.65 47.65 47.65 ... 47.95 47.95 47.95
    longitude  (pointid) float64 152B 15.99 16.19 16.39 16.03 ... 16.1 16.5 16.7
Data variables:
    2t         (t, pointid) float64 1MB 280.0 282.4 284.8 ... 276.6 272.6 273.9
    tp         (t, pointid) float64 1MB 0.0 0.0 0.0 ... 9.44e-06 3.725e-09 0.0

Plot both values over 2023

import matplotlib.pyplot as plt

# --- Plot 2t (2m Temperature) ---
fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True)

for point in ds_2023.pointid.values:
    axes[0].plot(
        ds_2023.t,
        ds_2023["2t"].sel(pointid=point) - 273.15,  # K to °C
        alpha=0.5, linewidth=0.7
    )
axes[0].set_title("2m Temperature – All Points (2023)")
axes[0].set_ylabel("Temperature (°C)")

# --- Plot tp (Total Precipitation) ---
for point in ds_2023.pointid.values:
    axes[1].plot(
        ds_2023.t,
        ds_2023["tp"].sel(pointid=point),
        alpha=0.5, linewidth=0.7
    )
axes[1].set_title("Total Precipitation – All Points (2023)")
axes[1].set_ylabel("Precipitation (m)")

plt.xlabel("Date")
plt.tight_layout()
plt.show()
../_images/df703f2a98071bbd5251b0189b8b8fe03e733da49ffe5111d857f16cb95de263.png