Get tiles from COGs

Demonstrates exploring HLS data at a specific time and place using the Raster API.
Author

Alexandra Kirk

Published

June 30, 2022

Run this notebook

You can launch this notbook using mybinder, by clicking the button below.

Binder

Approach

  1. Identify available dates within a bounding box, which is also an area of interest (AOI) in this example, for a given collection
  2. Register a dynamic tiler search for an AOI and specific date range for a given collection
  3. Explore different options for displaying multi-band Harmonized Landsat and Sentinel (HLS) assets with the Raster API.

About the Data

A small subset of HLS data has been ingested to the VEDA datastore to visually explore data using the Raster API, which is a VEDA instance of (pgstac-titiler). This limited subset includes a two granules for dates before and after Hurricane Maria in 2017 and Hurricane Ida in 2021.

Note about HLS datasets: The Sentinel and Landsat assets have been “harmonized” in the sense that these products have been generated to use the same spatial resolution and grid system. Thus these 2 HLS S30 and L30 productscan be used interchangeably in algorithms. However, the individual band assets are specific to each provider. This notebook focuses on displaying HLS data with a dynamic tiler so separate examples are provided for rendering the unique band assets of each collection.

Additional Resources

import json
import requests

from folium import Map, TileLayer

Parameters for investigating hurricane events with the dynamic tiler and custom band combinations

In this notebook we will focus on HLS S30 data for Hurricane Ida, but Hurricane Maria and L30 parameters are provided below for further exploration.

# Endpoints
STAC_API_URL = "https://staging-stac.delta-backend.com"
RASTER_API_URL = "https://staging-raster.delta-backend.com"

# Harmonized Sentinel collection id and configuration info
s30_collection_id = "hls-s30-002-ej-reprocessed"
s30_swir_assets = ["B12", "B8A", "B04"]
s30_vegetation_index_assets = ["B08", "B04"]
s30_vegetation_index_expression = "(B08_b1-B04_b1)/(B08_b1+B04_b1)"
s30_vegetation_index_rescaling = "0,1"
s30_vegetation_index_colormap = "rdylgn"

# Harmonized Landsat collection id and map configuration info
l30_collection_id = "hls-l30-002-ej-reprocessed"
l30_swir_assets = ["B07", "B05", "B04"]
l30_ndwi_expression = "(B03_b1-B05_b1)/(B03_b1+B05_b1)"
l30_ndwi_assets = ["B03", "B05"]
l30_ndwi_rescaling = "0,1"
l30_ndwi_colormap = "spectral"

# Search criteria for events in both HLS Events collections
maria_bbox = [-66.167596, 17.961538, -65.110098, 18.96772]
maria_temporal_range = ["2017-06-06T00:00:00Z", "2017-11-30T00:00:00Z"]

ida_bbox = [-90.932637, 29.705366, -89.766437, 30.71627]
ida_temporal_range = ["2021-07-01T00:00:00Z", "2021-10-28T00:00:00Z"]

First, search the STAC API to find the specific dates available within timeframe of interest (Hurricane Ida)

To focus on a specific point in time, we will restrict the temporal range when defining the item search in the example below.

collections_filter = {
    "op": "=",
    "args": [{"property": "collection"}, s30_collection_id],
}

spatial_filter = {"op": "s_intersects", "args": [{"property": "bbox"}, ida_bbox]}

temporal_filter = {
    "op": "t_intersects",
    "args": [{"property": "datetime"}, {"interval": ida_temporal_range}],
}

# Additional filters can be applied for other search criteria like <= maximum eo:cloud_cover in item properties
cloud_filter = {"op": "<=", "args": [{"property": "eo:cloud_cover"}, 80]}

search_body = {
    "filter-lang": "cql2-json",
    "limit": 100,
    "sortby": [{"direction": "asc", "field": "properties.datetime"}],
    "context": "on",  # add context for a summary of matched results
    "filter": {
        "op": "and",
        "args": [collections_filter, spatial_filter, temporal_filter, cloud_filter],
    },
}

# Note this search body can also be used for a stac item search
stac_items_response = requests.post(
    f"{STAC_API_URL}/search",
    json=search_body,
).json()

# Check how many items were matched in search
print("search context:", stac_items_response["context"])

# Iterate over search results to get an array of item datetimes
[item["properties"]["datetime"] for item in stac_items_response["features"]]
search context: {'limit': 100, 'matched': 14, 'returned': 14}
['2021-07-14T16:55:15.122720Z',
 '2021-07-24T16:55:15.112940Z',
 '2021-07-29T16:55:16.405890Z',
 '2021-08-08T16:55:15.798510Z',
 '2021-08-13T16:55:13.394950Z',
 '2021-08-23T16:55:11.785040Z',
 '2021-09-02T16:55:09.568600Z',
 '2021-09-07T16:55:13.430530Z',
 '2021-09-22T16:55:10.763010Z',
 '2021-09-27T16:55:17.027350Z',
 '2021-10-07T16:55:18.213640Z',
 '2021-10-12T16:55:14.209080Z',
 '2021-10-17T16:55:18.517600Z',
 '2021-10-22T16:55:14.670710Z']

Visualizing the data on a map

The VEDA backend is based on eoAPI, an application for searching and tiling earth observation STAC records. The application uses titiler-pgstac for dynamically mosaicing cloud optimized data from a registerd STAC API search.

To use the dynamic tiler, register a STAC item search and then use the registered search ID to dynamically mosaic the search results on the map.

Update the temporal range in search body and register that search with the Raster API

The registered search id can be reused for alternate map layer visualizations.

# Restricted date range
restricted_temporal_filter = {
    "op": "t_intersects",
    "args": [
        {"property": "datetime"},
        {"interval": ["2021-10-16T00:00:00Z", "2021-10-18T00:00:00Z"]},
    ],
}

# Specify cql2-json filter language in search body
search_body = {
    "filter-lang": "cql2-json",
    "filter": {
        "op": "and",
        "args": [collections_filter, spatial_filter, restricted_temporal_filter],
    },
}

mosaic_response = requests.post(
    f"{RASTER_API_URL}/mosaic/register",
    json=search_body,
).json()
print(json.dumps(mosaic_response, indent=2))
{
  "searchid": "7743bcb31bff7151aff7e5508785fce1",
  "links": [
    {
      "rel": "metadata",
      "title": "Mosaic metadata",
      "type": "application/json",
      "href": "https://staging-raster.delta-backend.com/mosaic/7743bcb31bff7151aff7e5508785fce1/info"
    },
    {
      "rel": "tilejson",
      "title": "Link for TileJSON",
      "type": "application/json",
      "href": "https://staging-raster.delta-backend.com/mosaic/7743bcb31bff7151aff7e5508785fce1/tilejson.json"
    },
    {
      "rel": "wmts",
      "title": "Link for WMTS",
      "type": "application/json",
      "href": "https://staging-raster.delta-backend.com/mosaic/7743bcb31bff7151aff7e5508785fce1/WMTSCapabilities.xml"
    }
  ]
}
# Get base url for tiler from the register mosaic request
tiles_href = next(
    link["href"] for link in mosaic_response["links"] if link["rel"] == "tilejson"
)

Configure map formatting parameters

See the raster-api/docs for more formatting options

Use the built-in SWIR post processing algorithm

Note in the example below the band assets for HLS S30 are selected. The equivalent SWIR band assets for L30 are provided at the top of this notebook.

# Add additional map formatting parameters to tiles url
tilejson_response = requests.get(
    tiles_href,
    params={
        # Info to add to the tilejson response
        "minzoom": 6,
        "maxzoom": 12,
        "algorith": "swir",
        "assets": s30_swir_assets,
    },
).json()
print(json.dumps(tilejson_response, indent=2))
{
  "tilejson": "2.2.0",
  "name": "7743bcb31bff7151aff7e5508785fce1",
  "version": "1.0.0",
  "scheme": "xyz",
  "tiles": [
    "https://staging-raster.delta-backend.com/mosaic/7743bcb31bff7151aff7e5508785fce1/tiles/WebMercatorQuad/{z}/{x}/{y}?algorith=swir&assets=B12&assets=B8A&assets=B04"
  ],
  "minzoom": 6,
  "maxzoom": 12,
  "bounds": [
    -180.0,
    -85.0511287798066,
    180.00000000000009,
    85.0511287798066
  ],
  "center": [
    4.263256414560601e-14,
    0.0,
    6
  ]
}

Display the data on a map

# Use bbox initial zoom and map
# Set up a map located w/in event bounds
zoom_start = 11
m = Map(
    tiles="OpenStreetMap",
    location=((ida_bbox[1] + ida_bbox[3]) / 2, (ida_bbox[0] + ida_bbox[2]) / 2),
    zoom_start=zoom_start,
)

# Add the formatted map layer
map_layer = TileLayer(
    tiles=tilejson_response["tiles"][0],
    attr="Mosaic",
)
map_layer.add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook

Format and render tiles using custom formatting

The titiler/raster-api supports user defined band combinations, band math expressions, rescaling, band index, resampling and more.

# Add additional map formatting parameters to tiles url
tilejson_response = requests.get(
    tiles_href,
    params={
        # Info to add to the tilejson response
        "minzoom": 6,
        "maxzoom": 12,
        "assets": s30_vegetation_index_assets,
        "expression": s30_vegetation_index_expression,
        "rescale": s30_vegetation_index_rescaling,
        "colormap_name": s30_vegetation_index_colormap,
    },
).json()
print(json.dumps(tilejson_response, indent=2))
{
  "tilejson": "2.2.0",
  "name": "7743bcb31bff7151aff7e5508785fce1",
  "version": "1.0.0",
  "scheme": "xyz",
  "tiles": [
    "https://staging-raster.delta-backend.com/mosaic/7743bcb31bff7151aff7e5508785fce1/tiles/WebMercatorQuad/{z}/{x}/{y}?assets=B08&assets=B04&expression=%28B08_b1-B04_b1%29%2F%28B08_b1%2BB04_b1%29&rescale=0%2C1&colormap_name=rdylgn"
  ],
  "minzoom": 6,
  "maxzoom": 12,
  "bounds": [
    -180.0,
    -85.0511287798066,
    180.00000000000009,
    85.0511287798066
  ],
  "center": [
    4.263256414560601e-14,
    0.0,
    6
  ]
}
# Use bbox initial zoom and map
# Set up a map located w/in event bounds
zoom_start = 11
m = Map(
    tiles="OpenStreetMap",
    location=((ida_bbox[1] + ida_bbox[3]) / 2, (ida_bbox[0] + ida_bbox[2]) / 2),
    zoom_start=zoom_start,
)

# Add the formatted map layer
map_layer = TileLayer(
    tiles=tilejson_response["tiles"][0],
    attr="Mosaic",
)
map_layer.add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook

L30 Hurricane Maria Example

collections_filter = {
    "op": "=", 
    "args" : [{ "property": "collection" }, l30_collection_id]
}

spatial_filter = {
    "op": "s_intersects",
    "args": [
        {"property": "bbox"}, maria_bbox
    ]
}

temporal_filter = {
    "op": "t_intersects",
    "args": [
        { "property": "datetime" },
        { "interval" : maria_temporal_range }
    ]
}

# Additional filters can be applied for other search criteria like <= maximum eo:cloud_cover in item properties
cloud_filter = {
    "op": "<=",
    "args": [
        {"property": "eo:cloud_cover"},
        80
    ]
}

# Specify cql2-json filter language in search body and add context for a summary of matched results
search_body = {
    "filter-lang": "cql2-json",
    "context": "on",
    "filter": {
        "op": "and",
        "args": [
            collections_filter,
            temporal_filter,
            cloud_filter
        ]
    }
}

# Note this search body can also be used for a stac item search 
stac_items_response = requests.post(
    f"{STAC_API_URL}/search",
    json=search_body,
).json()

# Check how many items were matched in searc
print("search context:", stac_items_response["context"])

# Iterate over search results to get an array of unique item datetimes
datetimes = []
features = stac_items_response["features"]
datetimes += [item["properties"]["datetime"] for item in features]
next_link = next((link for link in stac_items_response["links"] if link["rel"] == "next"), None)
while next_link:
    stac_items_response = requests.post(
        f"{STAC_API_URL}/search",
        json=next_link["body"],
    ).json()
    features = stac_items_response["features"]
    datetimes += [item["properties"]["datetime"] for item in features]
    next_link = next((link for link in stac_items_response["links"] if link["rel"] == "next"), False)

sorted(datetimes)
search context: {'limit': 10, 'matched': 9, 'returned': 9}
['2017-06-06T14:43:41.335694Z',
 '2017-06-22T14:43:47.156698Z',
 '2017-07-24T14:43:56.898518Z',
 '2017-08-09T14:44:03.584741Z',
 '2017-08-25T14:44:07.854507Z',
 '2017-09-26T14:44:14.813967Z',
 '2017-10-12T14:44:19.576858Z',
 '2017-11-13T14:44:17.834919Z',
 '2017-11-29T14:44:11.126689Z']
# Restricted date range 
restricted_temporal_filter = {
    "op": "t_intersects",
    "args": [
        { "property": "datetime" },
        { "interval" : [ "2017-10-11T00:00:00Z", "2017-10-13T00:00:00Z"] }
    ]
}

# Specify cql2-json filter language in search body
search_body = {
    "filter-lang": "cql2-json",
    "filter": {
        "op": "and",
        "args": [
            collections_filter,
            spatial_filter,
            restricted_temporal_filter
        ]
    }
}

mosaic_response = requests.post(
    f"{RASTER_API_URL}/mosaic/register",
    json=search_body,
).json()

# Set up format for Map API url
# Get base url for tiler from the register mosaic request
tiles_href = next(link["href"] for link in mosaic_response["links"] if link["rel"]=="tilejson")

# Add additional map formatting parameters to tiles url
tilejson_response = requests.get(
    tiles_href,
    params={
        # Info to add to the tilejson response
        "minzoom": 6,
        "maxzoom": 12,
        "assets": l30_ndwi_assets,
        "expression": l30_ndwi_expression,
        "rescale": l30_ndwi_rescaling,
        "colormap_name": "viridis"
    }
).json()
# Use bbox initial zoom and map
# Set up a map located w/in event bounds
zoom_start = 11
m = Map(
    tiles="OpenStreetMap",
    location=((maria_bbox[1] + maria_bbox[3]) / 2,(maria_bbox[0] + maria_bbox[2]) / 2),
    zoom_start=zoom_start
)

# Add the formatted map layer
map_layer = TileLayer(
    tiles=tilejson_response["tiles"][0],
    attr="Mosaic",  
)
map_layer.add_to(m)
m
Make this Notebook Trusted to load map: File -> Trust Notebook