Zarr Metadata With xstac - Ocean Hackathon

import xarray as xr
from rasterio.warp import transform_bounds
import rioxarray
from xstac import xarray_to_stac
import pandas as pd
import requests
from shapely.geometry import mapping, box

Open The Remote Zarr Store and Attach CRS Metadata¶

collection = 'sofresh-sea-surface-salinity'
zarr_href = 'https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/polar_cube_datasets/so-fresh/sofresh.zarr/'
ds = xr.open_zarr(zarr_href)
ds

ds = ds.rio.write_crs('ESRI:102020')
ds

Derive Spatial And Temporal Extents¶

This cell reads the projected x and y coordinate limits, transforms them to WGS84, then sets the final geographic bounding box used in STAC. It also converts the dataset time coordinate range into Python datetimes.

minx, maxx, miny, maxy = ds.x.min().values, ds.x.max().values, ds.y.min().values, ds.y.max().values

minx = float(ds.x.min())
maxx = float(ds.x.max())
miny = float(ds.y.min())
maxy = float(ds.y.max())

minx_4326, miny_4326, maxx_4326, maxy_4326 = transform_bounds(
    ds.rio.crs, "EPSG:4326", minx, miny, maxx, maxy, densify_pts=101
)

min_time, max_time = ds.time.values.ravel().min(), ds.time.values.ravel().max()


minx_4326, miny_4326, maxx_4326, maxy_4326 = (-180.0, -90.0, 180.0, -56.7587107166777)
start = pd.Timestamp(min_time).tz_localize("UTC").to_pydatetime()
end = pd.Timestamp(max_time).tz_localize("UTC").to_pydatetime()

Build The STAC Template And Let xstac Extract Metadata¶

The template supplies item-level fields such as ID, licence, temporal extent, geometry, and the Zarr data asset. xarray_to_stac then combines that template with metadata extracted from the xarray dataset. The final item.validate() call checks the generated STAC object with pystac validation.

raw_created = ds.attrs.get("product_created", ds.attrs.get("date_modified", pd.Timestamp.utcnow()))
created_utc = pd.to_datetime(raw_created, utc=True).isoformat().replace("+00:00", "Z")

start_utc = pd.to_datetime(start, utc=True).isoformat().replace("+00:00", "Z")
end_utc = pd.to_datetime(end, utc=True).isoformat().replace("+00:00", "Z")

item_id = f"{collection}-item"

template = {
    "id": item_id,
    "type": "Feature",
    "stac_version": "1.0.0",
    "properties": {
        "title": f"Sea surface salinity zarr store",
        "start_datetime": start_utc,
        "end_datetime": end_utc,
        "license": 'cc-by-4.0',
        "created": created_utc,
    },
    "bbox": [minx_4326, miny_4326, maxx_4326, maxy_4326],
    "geometry": mapping(box(minx_4326, miny_4326, maxx_4326, maxy_4326)),
    "assets": {
        "data": {
            "href": zarr_href,
            "type": "application/vnd+zarr",
            "roles": ["data"],
            "title": "Zarr Store",
        }
    },
}

item = xarray_to_stac(
    ds,
    template,
    temporal_dimension="time",
    x_dimension="x",
    y_dimension="y",
    validate=False
)

### bypass proj http block by m,anually downloading the schema
from pystac.validation import RegisteredValidator

url = 'https://proj.org/en/latest/schemas/v0.4/projjson.schema.json'
response = requests.get(url)
response.raise_for_status()
schema_data = response.json()

validator = RegisteredValidator.get_validator()
schema_url = "https://proj.org/schemas/v0.4/projjson.schema.json"
validator.schema_cache[schema_url] = schema_data


item.validate()

['https://schemas.stacspec.org/v1.1.0/item-spec/json-schema/item.json',
 'https://stac-extensions.github.io/datacube/v2.2.0/schema.json']