odc.stack.stac_load

Query Sentinel-2 data via odc.stac.stac_load into xarray.Dataset

from pystac_client import Client
import planetary_computer
from odc.stac import stac_load

from bounding_box import *

Important

Set modifier in client. Otherwise only the search will succeed, but data access will fail

catalog = Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

query = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[lon_min, lat_min, lon_max, lat_max],
    datetime="2024-07-01",
)

items = list(query.items())
print(f"Found: {len(items):d} datasets")

Important

Remember to set chunks, otherwise the data will be loaded eagerly!

chunk_size = 1024

ds = stac_load(
    items,
    bands=("B04", "B03", "B02", "B08"),
    resolution=10,
    chunks={"time": 1, "x": chunk_size, "y": chunk_size} # this prevents eager loading of data
)
ds

ds = ds.sel(x=slice(x_min, x_max), y=slice(y_max, y_min))

ds = ds.isel(time=[0]) # just take first timestep to reduce amount of data

# load data into memory if desired
ds.load()

ds["B04"].squeeze().plot.imshow()

… write to disk, if necessary …