Load Data from the eodc STAC catalogue into xarray

try:
    from odc import stac as odc_stac
except ImportError:
    %pip install odc-stac
    from odc import stac as odc_stac

try:
    from netCDF4 import Dataset
except ImportError:
    %pip install netCDF4
    from netCDF4 import Dataset

import xarray as xr
import matplotlib.pyplot as plt

from pystac_client import Client

import os
from datetime import datetime

import numpy as np

Access the catalogue and desired collection

First we find the data we want to load into xarray. This is described in the Data Discovery notebook.

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

time_range = "2024-03-01/2024-05-01"

# Bounding box of Austria
bbox_aut = [9.25, 46.31, 17.46, 49.18]

search = eodc_catalog.search(
    collections=["SENTINEL2_L1C"],
    bbox = bbox_aut,
    datetime=time_range
)

print("We found ", search.matched(), " items, that match our filter criteria.")

items = search.items()

We found  14  items, that match our filter criteria.

Specify Dataset parameters and load into xarray

crs = "EPSG:4312"
res = 0.00018 # 20 meters in degree 
bands = ["red", "nir"]

#This step can take time
print ("loading data started  : ", datetime.now())
ds = odc_stac.load(items,
                    crs=crs,
                    resolution=res,
                    bbox=bbox_aut,
                    bands = bands,
                       )
print ("loading data completed: ", datetime.now())

loading data started  :  2024-07-04 13:19:36.876273
loading data completed:  2024-07-04 13:49:39.271929

Save results as netCDF

ds.to_netcdf("./data/stac_xarray.nc")

Load netCDF

Open the saved NetCDF file.

nc_file = xr.open_dataset("./data/stac_xarray.nc")

red = nc_file['red']
nir = nc_file['nir']

def normalized_difference(x, y):
    nd = (x - y) / (x + y)
    return nd

ndvi = normalized_difference(nir, red)

Plot the NetCDF Data

Select data for a specific time step and create a plot.

time_data = ndvi.isel(time = 0)

p = plt.imshow(time_data, vmin = -1,  vmax = 1, cmap = "inferno")
plt.colorbar(p)
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("NetCDF Data NDVI")
plt.show()

To create different plots, you need to run the notebook yourself.