25  STAC specification

So far in our course we have accessed data in two ways: by downloading it directly from the data provider or by obtaining a URL from a data repository. This can be a convenient way to access targeted datasets, often usign GUIs (graphical user interfaces) for data discovery and filtering. However, relying on clicking and copy-pasting addresses and file names can make our workflows more error-prone and less reproducible. In particular, satellites around the world produce terabytes of new data daily and manually browsing through data repositories would it make hard to access this data. Moreover, we wouldn’t want to learn a new way to access data from every single big data provider. This is where STAC comes in.

The SpatioTemporal Asset Catalog (STAC) is an emerging open standard for geospatial data that aims to increase the interoperability of geospatial data, particularly satellite imagery. Many major data archives now follow the STAC specification.

In the next classes we’ll be working with the Microsoft’s Planetary Computer (MPC) STAC API. The MPC is both a geospatial coding environment and a STAC data catalog. In this lesson we will learn about the main components of a STAC catalog and how to search for data using the MPC’s STAC API.

25.1 Item, Collection, and Catalog

The STAC item (or just item) is the building block of a STAC. An item is a GeoJSON feature with additional fields that make it easier to find the item as we look for data across catalogs.

An item holds two types of information:

  1. Metadata: The metadata for a STAC item includes core identifying information (such as id, geometry, bounding box, and date), and additional properties (for example, place of collection).

  2. Assets: Assets are links to the actual data of the item (for example, links to the spectral bands of a satellite image.)

STAC items can be grouped into STAC collections. For example, while a single satellite scene (at a single time and location) would constitue an item, scenes across time and location from the same satellite can be orgnanized in a collection. Finally, multiple collections can be organized into a single STAC catalog.

For example, we’ll be accessing the Microsoft Planetary STAC catalog. Two of its collections are the National Agriculture Imagery Program (NAIP) colelction and the Copernicus Digital Elevation Model (DEM) colleciton. Each of these collections has multiple items, with item cotaining properties (metadata) and assets (links to the data).

25.2 API

To request data from a catalog following the STAC standard we use an Application Programming Interface (API). We can think of an API as an intermediary tasked with sending our request for data to the data catalog and getting the response from the catalog back to us. The following is a nice real-life analogy:

Geeks for geeks - What is an API?

The Python package to access APIs for STAC catalogs is pystac_client. Our goal in this lesson is to retrieve NAIP data from the MPC’s data catalog via its STAC API.

25.3 Catalog

First, load the necessary packages:

import numpy as np
import geopandas as gpd
import rioxarray as rioxr
import matplotlib.pyplot as plt

from shapely.geometry import Polygon

# used to access STAC catalogs
from pystac_client import Client
# used to sign items from the MPC STAC catalog
import planetary_computer

# ----- other libraries for nice ouputs
from IPython.display import Image

25.3.1 Access

We use the Client function from the pystac_client package to access the catalog:

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

The modifier parameter is needed to access the data in the MPC catalog.

25.3.2 Exploration

Let’s check out some of the catalog’s metadata:

# metadata from the catalog
print('Title:', catalog.title)
print('Description:', catalog.description)
Title: Microsoft Planetary Computer STAC API
Description: Searchable spatiotemporal metadata describing Earth science datasets hosted by the Microsoft Planetary Computer

We can access its collections by using the get_collections() method:

catalog.get_collections()
<generator object Client.get_collections at 0x112732570>

Notice the output of get_collections() is a generator. This is a special kind of lazy object in Python over which you can loop over like a list. Unlike a list, the items in a generator do not exist in memory until you explicitely iterate over them or convert them to a list. Let’s try getting the collections from the catalog again:

# get collections and print their names
collections = list(catalog.get_collections())

print('Number of collections:', len(collections))

print("Collections IDs (first 10):")
for i in range(10):
    print('-', collections[i].id)
Number of collections: 123
Collections IDs (first 10):
- daymet-annual-pr
- daymet-daily-hi
- 3dep-seamless
- 3dep-lidar-dsm
- fia
- sentinel-1-rtc
- gridmet
- daymet-annual-na
- daymet-monthly-na
- daymet-annual-hi

25.4 Collection

The NAIP catalog’s id is 'naip'. We can select a single collection for exploration using the get_child() method for the catalog and the collection id as the parameter:

naip_collection = catalog.get_child('naip')
naip_collection

25.6 Item

Let’s get the first item in the search:

# get first item in the catalog search
item = items[0]
type(item)
pystac.item.Item

Remember the STAC item is the core object in a STAC catalog. The item does not contain the data itself, but rather metadata and assets that contain links to access the actual data. Some of the metadata:

# print item id and properties
print('id:' , item.id)
item.properties
id: ca_m_3411935_sw_11_060_20200521
{'gsd': 0.6,
 'datetime': '2020-05-21T00:00:00Z',
 'naip:year': '2020',
 'proj:bbox': [246930.0, 3806808.0, 253260.0, 3814296.0],
 'proj:epsg': 26911,
 'naip:state': 'ca',
 'proj:shape': [12480, 10550],
 'proj:transform': [0.6, 0.0, 246930.0, 0.0, -0.6, 3814296.0, 0.0, 0.0, 1.0]}

Just as the item properties, the item assets are given in a dictionary, with each value being a pystac.asset Let’s check the assets in the item:

item.assets
{'image': <Asset href=https://naipeuwest.blob.core.windows.net/naip/v002/ca/2020/ca_060cm_2020/34119/m_3411935_sw_11_060_20200521.tif?st=2023-12-07T01%3A22%3A02Z&se=2023-12-08T02%3A07%3A02Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-12-08T01%3A19%3A05Z&ske=2023-12-15T01%3A19%3A05Z&sks=b&skv=2021-06-08&sig=KDgZ9RTP61bQDW47WsfPc/6czFOqQkR7FDSoCZigT8Y%3D>,
 'thumbnail': <Asset href=https://naipeuwest.blob.core.windows.net/naip/v002/ca/2020/ca_060cm_2020/34119/m_3411935_sw_11_060_20200521.200.jpg?st=2023-12-07T01%3A22%3A02Z&se=2023-12-08T02%3A07%3A02Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-12-08T01%3A19%3A05Z&ske=2023-12-15T01%3A19%3A05Z&sks=b&skv=2021-06-08&sig=KDgZ9RTP61bQDW47WsfPc/6czFOqQkR7FDSoCZigT8Y%3D>,
 'tilejson': <Asset href=https://planetarycomputer.microsoft.com/api/data/v1/item/tilejson.json?collection=naip&item=ca_m_3411935_sw_11_060_20200521&assets=image&asset_bidx=image%7C1%2C2%2C3&format=png>,
 'rendered_preview': <Asset href=https://planetarycomputer.microsoft.com/api/data/v1/item/preview.png?collection=naip&item=ca_m_3411935_sw_11_060_20200521&assets=image&asset_bidx=image%7C1%2C2%2C3&format=png>}
for key in item.assets.keys():
    print(key, '--', item.assets[key].title)
image -- RGBIR COG tile
thumbnail -- Thumbnail
tilejson -- TileJSON with default rendering
rendered_preview -- Rendered preview

Notice each asset has an href, which is a link to the data. For example, we can use the URL for the 'rendered_preview' asset to plot it:

# plot rendered preview
Image(url=item.assets['rendered_preview'].href, width=500)

25.7 Load data

The raster data in our current item is in the image asset. Again, we access this data via its URL. This time, we open it using rioxr.open_rasterio() directly:

sb = rioxr.open_rasterio(item.assets['image'].href)
sb
<xarray.DataArray (band: 4, y: 12480, x: 10550)>
[526656000 values with dtype=uint8]
Coordinates:
  * band         (band) int64 1 2 3 4
  * x            (x) float64 2.469e+05 2.469e+05 ... 2.533e+05 2.533e+05
  * y            (y) float64 3.814e+06 3.814e+06 ... 3.807e+06 3.807e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:             Area
    TIFFTAG_IMAGEDESCRIPTION:  OrthoVista
    TIFFTAG_RESOLUTIONUNIT:    1 (unitless)
    TIFFTAG_SOFTWARE:          Trimble Germany GmbH
    TIFFTAG_XRESOLUTION:       1
    TIFFTAG_YRESOLUTION:       1
    _FillValue:                0
    scale_factor:              1.0
    add_offset:                0.0

Notice this raster has four bands (red, green, blue, nir), so we cannot use the .plot.imshow() method directly (as this function only works when we have three bands). Thus we need select the bands we want to plot (RGB) before plotting:

# plot raster with correct ratio
size = 6 # height in in of plot height
aspect = sb.rio.width / sb.rio.height 
# select R,G,B bands and plot
sb.sel(band=[1,2,3]).plot.imshow(size=size, aspect=aspect)

25.8 Exercise

The 'cop-dem-glo-90' collection contains the Copernicus DEM at 90m resolution data (the data we previously used for the Grand Canyon).

  1. Reuse the bbox for Santa Barbara to look for items in this collection.
  2. Get the first item in the search and check its assets.
  3. Check the item’s rendered preview asset by clicking on it’s URL.
  4. Open the item’s data using rioxarray.

25.9 References

STAC Documentation:

Microsoft Planetary Computer Documentation - Reading Data from the STAC API