15 STAC specification

So far in our course we have obtained 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 graphical user interfaces (GUIs) 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 can it make difficult to access this data. Moreover, it can be inconvenient 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. 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.

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 it 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 MPC 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).

Application Programming Interface (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 diagram nicely explains what an API does using a real-life analogy of a restaurant:

Image Source: 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.

MPC Catalog

First, load the necessary packages:

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

from pystac_client import Client  # To access STAC catalogs

import planetary_computer  # To sign items from the MPC STAC catalog 

from IPython.display import Image  # To nicely display images

Access

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

# Access MPC 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.

Catalog Exploration

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

# Explore catalog metadata
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 0x1494fc590>

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. This can allow for more efficient memory allcoation. Once the generator is exhausted (i.e. iterated over completely), it cannot be reused unless it is recreated.

Let’s try getting the collections from the catalog again:

# Get collections and print their names
collections = list(catalog.get_collections())  # Turn generator into list

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

print("Collections IDs (first 10):")
for i in range(10):
    print('-', collections[i].id)
Number of collections: 125
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

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

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_20220513
{'gsd': 0.6,
 'datetime': '2022-05-13T16:00:00Z',
 'naip:year': '2022',
 'proj:bbox': [246930.0, 3806808.0, 253260.0, 3814296.0],
 'proj:epsg': 26911,
 'providers': [{'url': 'https://www.fsa.usda.gov/programs-and-services/aerial-photography/imagery-programs/naip-imagery/',
   'name': 'USDA Farm Service Agency',
   'roles': ['producer', 'licensor']}],
 'naip:state': 'ca',
 'proj:shape': [12480, 10550],
 'proj:centroid': {'lat': 34.40624, 'lon': -119.71877},
 '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/2022/ca_060cm_2022/34119/m_3411935_sw_11_060_20220513.tif?st=2024-11-25T22%3A06%3A49Z&se=2024-11-26T22%3A51%3A49Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-11-26T15%3A06%3A26Z&ske=2024-12-03T15%3A06%3A26Z&sks=b&skv=2024-05-04&sig=FDtjYiXqRfXkOsqqBYOJNvO4ozwvXKxaZEkS42czSWA%3D>,
 'thumbnail': <Asset href=https://naipeuwest.blob.core.windows.net/naip/v002/ca/2022/ca_060cm_2022/34119/m_3411935_sw_11_060_20220513.200.jpg?st=2024-11-25T22%3A06%3A49Z&se=2024-11-26T22%3A51%3A49Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-11-26T15%3A06%3A26Z&ske=2024-12-03T15%3A06%3A26Z&sks=b&skv=2024-05-04&sig=FDtjYiXqRfXkOsqqBYOJNvO4ozwvXKxaZEkS42czSWA%3D>,
 'tilejson': <Asset href=https://planetarycomputer.microsoft.com/api/data/v1/item/tilejson.json?collection=naip&item=ca_m_3411935_sw_11_060_20220513&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_20220513&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)

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
    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  
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)

Exercise

The 'cop-dem-glo-90' collection contains the Copernicus Digital Elevation Model (DEM) at 90m resolution data.

  1. Reuse the bbox for Santa Barbara to look for items in this collection.
  2. Get the first item in the search and examine its assets.
  3. Check the item’s rendered preview asset by clicking on it’s URL.
  4. Open and plot the item’s data using rioxarray.
  5. Obtain the maximum and minimum elevation on the scene as numbers.
  6. Print the maximum and minimum elevation rounded to two decimal points using f-strings.

References

STAC Documentation:

Microsoft Planetary Computer Documentation - Reading Data from the STAC API