Visualizing MAAP STAC Dataset with titiler-pgstac
Authors: Henry Rodman (DevSeed)
Date: September 5, 2024
Description: In this notebook, we visualize Cloud-Optimized GeoTIFFs (COGs) from the icesat2-boreal collection in MAAP’s STAC using the /collections
and /searches
endpoints in titiler-pgstac.maap-project.org.
Run This Notebook
To access and run this tutorial within MAAP’s Algorithm Development Environment (ADE), please refer to the “Getting started with the MAAP” section of our documentation.
Disclaimer: it is highly recommended to run a tutorial within MAAP’s ADE, which already includes the required software dependencies. Running the tutorial outside of the MAAP ADE may lead to errors. The pangeo
workspace in the MAAP ADE is already configured with all of the dependencies and access credentials.
Note About the Data
From the collection description: > Aboveground biomass density c.2020 gridded to 30m derived from ICESat-2, Harmonized Landsat-Sentinel2 and Copernicus DEM. Band 1 is the data band. Band 2 is the standard deviation.
These COGs are indexed in the MAAP STAC.
Additional Resources
Importing and Installing Packages
To be able to run this notebook you’ll need the following requirements:
folium
httpx
If the packages below are not installed already, uncomment the following cell:
[1]:
# %pip install folium
# %pip install httpx
[2]:
import httpx
from folium import GeoJson, LayerControl, Map, TileLayer
from pprint import pprint
from pystac_client import Client
from shapely.geometry import box, mapping
This demo uses the MAAP titiler-pgstac deployment to render tiles for a collection (icesat2-boreal
) that is indexed in the MAAP STAC. We will make use of some pre-defined render
parameters stored in the `render
extension <https://github.com/stac-extensions/render>`__ metadata for the collection when we generate visualizations.
[3]:
TITILER_PGSTAC_ENDPOINT = "https://titiler-pgstac.maap-project.org" # MAAP titiler-pgstac endpoint
COLLECTION_ID = "icesat2-boreal"
RENDERING = "agb" # render parameters for aboveground biomass
[4]:
# get some information about the collection including render parameters and spatial extent
collection_info = httpx.get(
f"{TITILER_PGSTAC_ENDPOINT}/collections/{COLLECTION_ID}/info",
timeout=10,
).json()
pprint(collection_info)
{'links': [{'href': 'https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/info',
'rel': 'self',
'title': 'Mosaic metadata'},
{'href': 'https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/tilejson.json',
'rel': 'tilejson',
'templated': True,
'title': 'TileJSON link (Template URL).'},
{'href': 'https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/tilejson.json?bidx=1&nodata=nan&colormap_name=gist_earth_r&rescale=0%2C400&assets=tif',
'rel': 'tilejson',
'templated': True,
'title': 'TileJSON link for `agb` layer (Template URL).'},
{'href': 'https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/map',
'rel': 'map',
'templated': True,
'title': 'Map viewer link (Template URL).'},
{'href': 'https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/map?bidx=1&nodata=nan&colormap_name=gist_earth_r&rescale=0%2C400&assets=tif',
'rel': 'map',
'templated': True,
'title': 'Map viewer link for `agb` layer (Template URL).'},
{'href': 'https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/WMTSCapabilities.xml',
'rel': 'wmts',
'templated': True,
'title': 'WMTS link (Template URL)'}],
'search': {'_where': "collection = ANY ('{icesat2-boreal}') ",
'hash': '2b1c56237655c787dc5598c8592f5010',
'lastused': '2024-09-06T20:29:27.945205Z',
'metadata': {'assets': ['csv', 'tif'],
'bounds': [-180.0, 51.6, 180.0, 78.0],
'defaults': {'agb': {'assets': ['tif'],
'bidx': [1],
'colormap_name': 'gist_earth_r',
'nodata': 'nan',
'rescale': [[0, 400]]}},
'name': "Mosaic for 'icesat2-boreal' Collection",
'type': 'mosaic'},
'orderby': 'datetime DESC, id DESC',
'search': {'collections': ['icesat2-boreal']},
'usecount': 462}}
Tip: You can also get the render parameters from the /collections/{collection_id} STAC API endpoint but for simplicity we are using the titiler-pgstac endpoints for everything in this example.
Render tiles for the whole collection
The `/collections
endpoints <https://titiler-pgstac.maap-project.org/api.html#/STAC%20Collection>`__ in a titiler-pgstac
deployment can be used to render tiles for a mosaic of all items in a collection. This is particularly convenient for collections where all areas are covered by only a single item because you can get everything you need for rendering snazzy tiles with a single GET
request. The catch is that if your collection contains data for multiple time points, for example, the
tiles will load the most recent ones by default and stop adding data once the tile is fully covered by items (in descending order by time).
The render
parameters are provided in the collection metadata and can be accessed in the response from the /collections/{collection_id}/info
endpoint:
[5]:
# load the render parameters from the collection metadata
agb_render_params = collection_info["search"]["metadata"]["defaults"][RENDERING].copy()
# add some zoom limits to avoid overloading the tile server
agb_render_params["minzoom"] = 7
agb_render_params["maxzoom"] = 14
pprint(agb_render_params)
{'assets': ['tif'],
'bidx': [1],
'colormap_name': 'gist_earth_r',
'maxzoom': 14,
'minzoom': 7,
'nodata': 'nan',
'rescale': [[0, 400]]}
These parameters can be passed along to any titiler-pgstac
endpoint to make use of pre-defined visualization settings!
[6]:
collection_tilejson = httpx.get(
f"{TITILER_PGSTAC_ENDPOINT}/collections/{COLLECTION_ID}/WebMercatorQuad/tilejson.json",
params=agb_render_params,
).json()
pprint(collection_tilejson)
{'bounds': [-180.0, 51.6, 180.0, 78.0],
'center': [0.0, 64.8, 7],
'maxzoom': 14,
'minzoom': 7,
'name': "Mosaic for 'icesat2-boreal' Collection",
'scheme': 'xyz',
'tilejson': '2.2.0',
'tiles': ['https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/tiles/WebMercatorQuad/{z}/{x}/{y}?bidx=1&assets=tif&nodata=nan&rescale=%5B0%2C+400%5D&colormap_name=gist_earth_r'],
'version': '1.0.0'}
The contents of the collection_tilejson
dictionary can be passed to the arguments in TileLayer
to control things like min/max zoom levels for the tiles.
[7]:
# create an XYZ tile layer for the leaflet map
tiles = TileLayer(
name=f"{COLLECTION_ID} - {RENDERING}",
tiles=collection_tilejson["tiles"][0],
min_zoom=collection_tilejson["minzoom"],
max_zoom=collection_tilejson["maxzoom"],
opacity=1,
attr="NASA",
overlay=True,
control=True,
)
# load the collection bounding box so we can plot it on the map
collection_bbox = collection_info["search"]["metadata"]["bounds"]
geojson = {"type": "Feature", "geometry": mapping(box(*collection_bbox)), "properties": {}}
geojson_layer = GeoJson(
data=geojson,
style_function=lambda x: {
"opacity": 1,
"dashArray": "1",
"fillOpacity": 0,
"weight": 2,
"color": "orange",
},
name=f"{COLLECTION_ID} bounds",
overlay=True,
control=True,
show=True,
)
m = Map(
tiles="OpenStreetMap",
location=((collection_bbox[1] + collection_bbox[3]) / 2, (collection_bbox[0] + collection_bbox[2]) / 2),
zoom_start=2,
)
tiles.add_to(m)
geojson_layer.add_to(m)
LayerControl(collapsed=False).add_to(m)
m
[7]:
Render tiles for a STAC search
You can use the /searches
endpoints to register a STAC search with titiler-pgstac
and view tiles from a mosaic formed by the items that are returned by the search. This is useful for limiting tiles to a particular temporal or spatial extent or possibly to items with specific properties.
This can be done with a POST
request to the `/searches/register
endpoint <https://titiler-pgstac.maap-project.org/api.html#/STAC%20Search/register_search_searches_register_post>`__ with the STAC search parameters included as a json request body. You can use any search
parameters that are implemented in the pgstac
database backing the titiler-pgstac
deployment.
In this example, we will register a STAC search for items that intersect the bounding box of Alaska:
[8]:
ak_bbox = (-180, 50, -128, 73)
ak_search_request = httpx.post(
f"{TITILER_PGSTAC_ENDPOINT}/searches/register",
json={
"collections": [COLLECTION_ID],
"bbox": ak_bbox,
},
).json()
pprint(ak_search_request)
{'id': '058af6b76d8fc4e092818a0d1f7bffa8',
'links': [{'href': 'https://titiler-pgstac.maap-project.org/searches/058af6b76d8fc4e092818a0d1f7bffa8/info',
'rel': 'metadata',
'title': 'Mosaic metadata'},
{'href': 'https://titiler-pgstac.maap-project.org/searches/058af6b76d8fc4e092818a0d1f7bffa8/{tileMatrixSetId}/tilejson.json',
'rel': 'tilejson',
'templated': True,
'title': 'Link for TileJSON (Template URL)'},
{'href': 'https://titiler-pgstac.maap-project.org/searches/058af6b76d8fc4e092818a0d1f7bffa8/{tileMatrixSetId}/map',
'rel': 'map',
'templated': True,
'title': 'Link for Map viewer (Template URL)'},
{'href': 'https://titiler-pgstac.maap-project.org/searches/058af6b76d8fc4e092818a0d1f7bffa8/{tileMatrixSetId}/WMTSCapabilities.xml',
'rel': 'wmts',
'templated': True,
'title': 'Link for WMTS (Template URL)'}]}
The /searches
POST request will return a search id
that represents a unique key for the search parameters that you provided. The id
can be used in the /searches/{search_id}
endpoints! However, the render
parameters from the collection metadata are not be associated with any search
that is registered in this way so we need to re-use the agb_render_params
that we pulled from the /collections/{collection_id}/info
endpoint.
[9]:
search_id = ak_search_request["id"]
ak_tilejson = httpx.get(
f"{TITILER_PGSTAC_ENDPOINT}/searches/{search_id}/WebMercatorQuad/tilejson.json",
params=agb_render_params,
).json()
pprint(ak_tilejson)
{'bounds': [-180.0, 50.0, -128.0, 73.0],
'center': [-154.0, 61.5, 7],
'maxzoom': 14,
'minzoom': 7,
'name': '058af6b76d8fc4e092818a0d1f7bffa8',
'scheme': 'xyz',
'tilejson': '2.2.0',
'tiles': ['https://titiler-pgstac.maap-project.org/searches/058af6b76d8fc4e092818a0d1f7bffa8/tiles/WebMercatorQuad/{z}/{x}/{y}?bidx=1&assets=tif&nodata=nan&rescale=%5B0%2C+400%5D&colormap_name=gist_earth_r'],
'version': '1.0.0'}
[10]:
tiles = TileLayer(
name=f"{COLLECTION_ID} - {RENDERING}",
tiles=ak_tilejson["tiles"][0],
min_zoom=ak_tilejson["minzoom"],
max_zoom=ak_tilejson["maxzoom"],
opacity=1,
attr="NASA",
overlay=True,
control=True,
)
geojson = {"type": "Feature", "geometry": mapping(box(*ak_bbox)), "properties": {}}
zoom_start = 5
geojson_layer = GeoJson(
data=geojson,
style_function=lambda x: {
"opacity": 1,
"dashArray": "1",
"fillOpacity": 0,
"weight": 1,
},
name="STAC search bounds",
overlay=True,
control=True,
show=True,
)
m = Map(
tiles="OpenStreetMap",
location=((ak_bbox[1] + ak_bbox[3]) / 2, (ak_bbox[0] + ak_bbox[2]) / 2),
zoom_start=zoom_start,
)
geojson_layer.add_to(m)
tiles.add_to(m)
LayerControl(collapsed=False).add_to(m)
m
[10]:
[ ]: