How to read data from STAC#

This notebook shows how to read the data in from a STAC asset using xarray and a little hidden helper library called xpystac.

tl;dr#

For any PySTAC asset that can be represented as an xarray dataset you can read the data using the following command:

xr.open_dataset(asset)

If you want to load multiple assets from the same item and/or many items at once use:

odc.stac.load([items])

Dependencies#

There are lots of optional dependencies depending on where and how the data you are interested in are stored. Here are some of the libraries that you will probably need:

  • dask - to delay data loading until access

  • pystac - STAC object structures

  • xarray, rioxarray - data structures

  • xpystac, odc-stac - helper for loading pystac into xarray objects

Despite all these install instructions, the import block is very straightforward

[1]:

import odc.stac
import planetary_computer
import rioxarray
import xarray as xr

import pystac

Examples#

Here are a few examples of the different types of objects that you can open in xarray.

COGs#

Read all the data from the COGs referenced by the assets on an item.

[2]:

landsat_item = pystac.Item.from_file(
    "https://planetarycomputer.microsoft.com/api/stac/v1/collections/landsat-c2-l2/items/LC09_L2SP_088084_20230408_02_T2",
)

ds = odc.stac.load(
    [landsat_item], chunks={"x": 1024, "y": 1024}, patch_url=planetary_computer.sign
)
ds

[2]:

<xarray.Dataset> Size: 2GB
Dimensions:      (y: 7801, x: 7761, time: 1)
Coordinates:
  * y            (y) float64 62kB -3.713e+06 -3.713e+06 ... -3.947e+06
  * x            (x) float64 62kB 3.774e+05 3.774e+05 ... 6.102e+05 6.102e+05
    spatial_ref  int32 4B 32656
  * time         (time) datetime64[ns] 8B 2023-04-08T23:37:51.630731
Data variables: (12/19)
    qa           (time, y, x) int16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    red          (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    blue         (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    drad         (time, y, x) int16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    emis         (time, y, x) int16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    emsd         (time, y, x) int16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    ...           ...
    swir16       (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    swir22       (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    coastal      (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    qa_pixel     (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    qa_radsat    (time, y, x) uint16 121MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>
    qa_aerosol   (time, y, x) uint8 61MB dask.array<chunksize=(1, 1024, 1024), meta=np.ndarray>

Let’s prove that we really can access the data within the COGs.

[3]:

%%time
ds.blue.mean().compute()

CPU times: user 2.68 s, sys: 459 ms, total: 3.14 s
Wall time: 5.31 s

[3]:

<xarray.DataArray 'blue' ()> Size: 8B
array(13940.01841915)
Coordinates:
    spatial_ref  int32 4B 32656

For more control over the individual COGs referenced by the assets, you can grab the href off the asset, sign it and then use it directly.

[4]:

landsat_asset = landsat_item.assets["blue"]
landsat_asset_href = planetary_computer.sign(landsat_asset).href

da = rioxarray.open_rasterio(landsat_asset_href, overview_level=5)
da

[4]:

<xarray.DataArray (band: 1, y: 122, x: 122)> Size: 30kB
[14884 values with dtype=uint16]
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 976B 3.783e+05 3.802e+05 ... 6.074e+05 6.093e+05
  * y            (y) float64 976B -3.714e+06 -3.716e+06 ... -3.946e+06
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Point
    _FillValue:     0
    scale_factor:   1.0
    add_offset:     0.0
[5]:

%%time
da.mean()

CPU times: user 3.44 ms, sys: 0 ns, total: 3.44 ms
Wall time: 186 ms

[5]:

<xarray.DataArray ()> Size: 8B
array(14472.21996775)
Coordinates:
    spatial_ref  int64 8B 0

Zarr#

Read from an asset that references data stored in zarr

[6]:

daymet_collection = pystac.Collection.from_file(
    "https://planetarycomputer.microsoft.com/api/stac/v1/collections/daymet-daily-hi"
)
daymet_asset = daymet_collection.assets["zarr-abfs"]
daymet_asset_signed = planetary_computer.sign(daymet_asset)

ds = xr.open_dataset(daymet_asset_signed, chunks={})
ds

[6]:

<xarray.Dataset> Size: 69GB
Dimensions:                  (time: 14965, y: 584, x: 284, nv: 2)
Coordinates:
    lat                      (y, x) float32 663kB dask.array<chunksize=(584, 284), meta=np.ndarray>
    lon                      (y, x) float32 663kB dask.array<chunksize=(584, 284), meta=np.ndarray>
  * time                     (time) datetime64[ns] 120kB 1980-01-01T12:00:00 ...
  * x                        (x) float32 1kB -5.802e+06 ... -5.519e+06
  * y                        (y) float32 2kB -3.9e+04 -4e+04 ... -6.22e+05
Dimensions without coordinates: nv
Data variables:
    dayl                     (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    lambert_conformal_conic  int16 2B ...
    prcp                     (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    srad                     (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    swe                      (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    time_bnds                (time, nv) datetime64[ns] 239kB dask.array<chunksize=(365, 2), meta=np.ndarray>
    tmax                     (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    tmin                     (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    vp                       (time, y, x) float32 10GB dask.array<chunksize=(365, 584, 284), meta=np.ndarray>
    yearday                  (time) int16 30kB dask.array<chunksize=(365,), meta=np.ndarray>
Attributes:
    Conventions:       CF-1.6
    Version_data:      Daymet Data Version 4.0
    Version_software:  Daymet Software Version 4.0
    citation:          Please see http://daymet.ornl.gov/ for current Daymet ...
    references:        Please see http://daymet.ornl.gov/ for current informa...
    source:            Daymet Software Version 4.0
    start_year:        1980

Reference file#

If the collection has a reference file we can use that.

This will not work for kerchunk>=0.2.8 and xpystac<=0.2.0. Track the xpystac ticket here: xpystac #48

[7]:

cmip6_collection = pystac.Collection.from_file(
    "https://planetarycomputer.microsoft.com/api/stac/v1/collections/nasa-nex-gddp-cmip6"
)
cmip6_asset = cmip6_collection.assets["ACCESS-CM2.historical"]
cmip6_asset_signed = planetary_computer.sign(cmip6_asset)

ds = xr.open_dataset(cmip6_asset_signed, chunks={}, patch_url=planetary_computer.sign)
ds

[7]:

<xarray.Dataset> Size: 738GB
Dimensions:  (time: 23741, lat: 600, lon: 1440)
Coordinates:
  * lat      (lat) float64 5kB -59.88 -59.62 -59.38 -59.12 ... 89.38 89.62 89.88
  * lon      (lon) float64 12kB 0.125 0.375 0.625 0.875 ... 359.4 359.6 359.9
  * time     (time) datetime64[us] 190kB 1950-01-01T12:00:00 ... 2014-12-31T1...
Data variables:
    hurs     (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    huss     (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    pr       (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    rlds     (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    rsds     (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    sfcWind  (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    tas      (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    tasmax   (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
    tasmin   (time, lat, lon) float32 82GB dask.array<chunksize=(1, 600, 1440), meta=np.ndarray>
Attributes: (12/22)
    Conventions:           CF-1.7
    activity:              NEX-GDDP-CMIP6
    cmip6_institution_id:  CSIRO-ARCCSS
    cmip6_license:         CC-BY-SA 4.0
    cmip6_source_id:       ACCESS-CM2
    contact:               Dr. Rama Nemani: rama.nemani@nasa.gov, Dr. Bridget...
    ...                    ...
    scenario:              historical
    source:                BCSD
    title:                 ACCESS-CM2, r1i1p1f1, historical, global downscale...
    tracking_id:           16d27564-470f-41ea-8077-f4cc3efa5bfe
    variant_label:         r1i1p1f1
    version:               1.0