Working With STAC¶

STAC: SpatioTemporal Asset Catalog¶

The SpatioTemporal Asset Catalog (STAC) specification aims to standardize the way geospatial assets are exposed online and queried. A 'spatiotemporal asset' is any file that represents information about the earth captured in a certain space and time. The initial focus is primarily remotely-sensed imagery (from satellites, but also planes, drones, balloons, etc), but the core is designed to be extensible to SAR, full motion video, point clouds, hyperspectral, LiDAR and derived data like NDVI, Digital Elevation Models, mosaics, etc.

Ref: https://github.com/radiantearth/stac-spechttps://github.com/radiantearth/stac-spec

Using STAC makes data indexation and discovery really easy. In addition to the Collection/Item/Asset (data) specifications, data providers are also encouraged to follow a STAC API specification: https://github.com/radiantearth/stac-api-spec

The API is compliant with the OGC API - Features standard (formerly known as OGC Web Feature Service 3), in that it defines many of the endpoints that STAC uses. A STAC API should be compatible and usable with any OGC API - Features clients. The STAC API can be thought of as a specialized Features API to search STAC Catalogs, where the features returned are STAC Items, that have common properties, links to their assets and geometries that represent the footprints of the geospatial assets.

Requirements¶

To be able to run this notebook you'll need the following requirements:

• folium
• httpx
• rasterio

!pip install folium httpx rasterio

In [ ]:

# Uncomment this line if you need to install the dependencies
# !pip install folium requests rasterio

# Uncomment this line if you need to install the dependencies # !pip install folium requests rasterio

In [ ]:

import httpx

from rasterio.features import bounds as featureBounds

from folium import Map, TileLayer, GeoJson

%pylab inline

import httpx from rasterio.features import bounds as featureBounds from folium import Map, TileLayer, GeoJson %pylab inline

In [ ]:

titiler_endpoint = "https://titiler.xyz"  # Developmentseed Demo endpoint. Please be kind.
stac_item = "https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2A_30TVT_20221112_0_L2A"

titiler_endpoint = "https://titiler.xyz" # Developmentseed Demo endpoint. Please be kind. stac_item = "https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2A_30TVT_20221112_0_L2A"

In [ ]:

item = httpx.get(stac_item).json()
print(item)

item = httpx.get(stac_item).json() print(item)

In [ ]:

for it, asset in item["assets"].items():
    print("Name:", it, "| Format:", asset["type"])

for it, asset in item["assets"].items(): print("Name:", it, "| Format:", asset["type"])

In [ ]:

bounds = featureBounds(item)

m = Map(
    tiles="OpenStreetMap",
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=8
)

geo_json = GeoJson(data=item)
geo_json.add_to(m)
m

bounds = featureBounds(item) m = Map( tiles="OpenStreetMap", location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=8 ) geo_json = GeoJson(data=item) geo_json.add_to(m) m

In [ ]:

# Get Tile URL
r = httpx.get(
    f"{titiler_endpoint}/stac/info",
    params = (
        ("url", stac_item),
        # Get info for multiple assets
        ("assets","visual"), ("assets","red"), ("assets","blue"), ("assets","green"),
    )
).json()
print(r)

# Get Tile URL r = httpx.get( f"{titiler_endpoint}/stac/info", params = ( ("url", stac_item), # Get info for multiple assets ("assets","visual"), ("assets","red"), ("assets","blue"), ("assets","green"), ) ).json() print(r)

Display one asset¶

In [ ]:

r = httpx.get(
    f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json",
    params = {
        "url": stac_item,
        "assets": "visual",
        "minzoom": 8,  # By default titiler will use 0
        "maxzoom": 14, # By default titiler will use 24
    }
).json()

m = Map(
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=10
)

tiles = TileLayer(
    tiles=r["tiles"][0],
    min_zoom=r["minzoom"],
    max_zoom=r["maxzoom"],
    opacity=1,
    attr="ESA"
)
tiles.add_to(m)
m

r = httpx.get( f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json", params = { "url": stac_item, "assets": "visual", "minzoom": 8, # By default titiler will use 0 "maxzoom": 14, # By default titiler will use 24 } ).json() m = Map( location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=10 ) tiles = TileLayer( tiles=r["tiles"][0], min_zoom=r["minzoom"], max_zoom=r["maxzoom"], opacity=1, attr="ESA" ) tiles.add_to(m) m

Select Indexes for assets¶

In [ ]:

# Get Tile URL
r = httpx.get(
    f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json",
    params = {
        "url": stac_item,
        "assets": "visual",
        "asset_bidx": "visual|3,1,2",
        "minzoom": 8,  # By default titiler will use 0
        "maxzoom": 14, # By default titiler will use 24
    }
).json()

m = Map(
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=12
)

tiles = TileLayer(
    tiles=r["tiles"][0],
    min_zoom=r["minzoom"],
    max_zoom=r["maxzoom"],
    opacity=1,
    attr="ESA"
)
tiles.add_to(m)
m

# Get Tile URL r = httpx.get( f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json", params = { "url": stac_item, "assets": "visual", "asset_bidx": "visual|3,1,2", "minzoom": 8, # By default titiler will use 0 "maxzoom": 14, # By default titiler will use 24 } ).json() m = Map( location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=12 ) tiles = TileLayer( tiles=r["tiles"][0], min_zoom=r["minzoom"], max_zoom=r["maxzoom"], opacity=1, attr="ESA" ) tiles.add_to(m) m

In [ ]:

# Get Tile URL
r = httpx.get(
    f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json",
    params = (
        ("url", stac_item),
        ("assets", "red"),
        ("assets", "green"),
        ("assets", "blue"),
        # Most of the Sentinel L2A Assets have only one band
        # So we don't have to pass the bidx
        # ("assets_bidx", "red|1"),
        # ("assets_bidx", "green|1"),
        # ("assets_bidx", "blue|"),
        ("minzoom", 8),
        ("maxzoom", 14),
        ("rescale", "0,2000"),
    )
).json()

m = Map(
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=11
)

tiles = TileLayer(
    tiles=r["tiles"][0],
    min_zoom=r["minzoom"],
    max_zoom=r["maxzoom"],
    opacity=1,
    attr="ESA"
)
tiles.add_to(m)
m

# Get Tile URL r = httpx.get( f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json", params = ( ("url", stac_item), ("assets", "red"), ("assets", "green"), ("assets", "blue"), # Most of the Sentinel L2A Assets have only one band # So we don't have to pass the bidx # ("assets_bidx", "red|1"), # ("assets_bidx", "green|1"), # ("assets_bidx", "blue|"), ("minzoom", 8), ("maxzoom", 14), ("rescale", "0,2000"), ) ).json() m = Map( location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=11 ) tiles = TileLayer( tiles=r["tiles"][0], min_zoom=r["minzoom"], max_zoom=r["maxzoom"], opacity=1, attr="ESA" ) tiles.add_to(m) m

Use an expression to calculate a band index (NDVI) based on information contained in multiple assets.

In [ ]:

# Get Tile URL
r = httpx.get(
    f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json",
    params = (
        ("url", stac_item),
        ("expression", "(nir-red)/(nir+red)"),  # NDVI
        # We need to tell rio-tiler that each asset is a Band
        # (so it will select the first band within each asset automatically)
        ("asset_as_band", True),
        ("rescale", "-1,1"),
        ("minzoom", 8),
        ("maxzoom", 14),
        ("colormap_name", "viridis"),
    )
).json()

m = Map(
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=10
)

tiles = TileLayer(
    tiles=r["tiles"][0],
    min_zoom=r["minzoom"],
    max_zoom=r["maxzoom"],
    opacity=1,
    attr="ESA"
)
tiles.add_to(m)
m

# Get Tile URL r = httpx.get( f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json", params = ( ("url", stac_item), ("expression", "(nir-red)/(nir+red)"), # NDVI # We need to tell rio-tiler that each asset is a Band # (so it will select the first band within each asset automatically) ("asset_as_band", True), ("rescale", "-1,1"), ("minzoom", 8), ("maxzoom", 14), ("colormap_name", "viridis"), ) ).json() m = Map( location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=10 ) tiles = TileLayer( tiles=r["tiles"][0], min_zoom=r["minzoom"], max_zoom=r["maxzoom"], opacity=1, attr="ESA" ) tiles.add_to(m) m

If you don't use the asset_as_band=True option, you need to append the band to the asset name within the expression. For example, nir becomes nir_b1.

In [ ]:

# Get Tile URL
r = httpx.get(
    f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json",
    params = (
        ("url", stac_item),
        ("expression", "(nir_b1-red_b1)/(nir_b1+red_b1)"),  # NDVI
        ("rescale", "-1,1"),
        ("minzoom", 8),
        ("maxzoom", 14),
        ("colormap_name", "viridis"),
    )
).json()

m = Map(
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=10
)

tiles = TileLayer(
    tiles=r["tiles"][0],
    min_zoom=r["minzoom"],
    max_zoom=r["maxzoom"],
    opacity=1,
    attr="ESA"
)
tiles.add_to(m)
m

# Get Tile URL r = httpx.get( f"{titiler_endpoint}/stac/WebMercatorQuad/tilejson.json", params = ( ("url", stac_item), ("expression", "(nir_b1-red_b1)/(nir_b1+red_b1)"), # NDVI ("rescale", "-1,1"), ("minzoom", 8), ("maxzoom", 14), ("colormap_name", "viridis"), ) ).json() m = Map( location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=10 ) tiles = TileLayer( tiles=r["tiles"][0], min_zoom=r["minzoom"], max_zoom=r["maxzoom"], opacity=1, attr="ESA" ) tiles.add_to(m) m

In [ ]: