Visualizing MAAP STAC Dataset with titiler-pgstac in R

Authors: Harshini Girish (UAH), Sheyenne Kirkland (UAH), Alex Mandel (Development Seed), Henry Rodman (Development Seed)

Date: July 16, 2025

Description: n 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 packages specific to MAAP, such as maap-py. Running the tutorial outside of the MAAP ADE may lead to errors. Users should work within the “R/Python” workspace.

Additional Resources

Visualizing STAC Items: Official MAAP tutorial demonstrating how to visualize spatial footprints from STAC collections using different functions inside the ADE.
Titiler-PgSTAC API Reference: Full API reference for interacting with the Titiler-PgSTAC service, which enables dynamic tiling and rendering of geospatial assets from MAAP’s STAC catalog.
Visualizing with Titiler-PgSTAC: MAAP tutorial on rendering and visualizing COGs using the Titiler-PgSTAC tiling service, including usage examples and parameter explanations.

Install/Import Libraries

Let’s install and load the packages necessary for this tutorial.

[1]:

library(httr2)
library(jsonlite)
library(leaflet)
library(rstac)
library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

Fetch Collection Info from Endpoint

This code queries the /info endpoint for the icesat2-boreal collection using httr2. It retrieves metadata like bounds, mosaic type, and rendering defaults for use in tile visualization. The parsed JSON includes search settings, usage stats, and the collection name.

This information is used to configure how the collection is rendered as a map layer.

[2]:

titiler_endpoint <- "https://titiler-pgstac.maap-project.org"
collection_id <- "icesat2-boreal"
info_url <- paste0(titiler_endpoint, "/collections/", collection_id, "/info")

[3]:

response <- request(info_url) |> req_timeout(10) |> req_perform()
collection_info <- fromJSON(resp_body_string(response), simplifyVector = TRUE)
print(collection_info)

$search
$search$hash
[1] "e1458a47fc19d7e34c06f3d8469fd4f6"

$search$search
$search$search$collections
[1] "icesat2-boreal"

$search$search$`filter-lang`
[1] "cql2-json"


$search$lastused
[1] "2025-12-18T20:40:25.515213Z"

$search$usecount
[1] 175

$search$metadata
$search$metadata$type
[1] "mosaic"

$search$metadata$bounds
[1] -180.0   51.6  180.0   78.0

$search$metadata$name
[1] "Mosaic for 'icesat2-boreal' Collection"

$search$metadata$assets
[1] "csv" "tif"

$search$metadata$defaults
$search$metadata$defaults$agb
$search$metadata$defaults$agb$bidx
[1] 1

$search$metadata$defaults$agb$assets
[1] "tif"

$search$metadata$defaults$agb$nodata
[1] "nan"

$search$metadata$defaults$agb$rescale
     [,1] [,2]
[1,]    0  400

$search$metadata$defaults$agb$colormap_name
[1] "gist_earth_r"





$links
                                                                                                                                                                        href
1                                                                                                    https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/info
2                                                                         https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/tilejson.json
3 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
4                                                                              https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/map.html
5      https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/map.html?bidx=1&nodata=nan&colormap_name=gist_earth_r&rescale=0%2C400&assets=tif
6                                                                  https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/WMTSCapabilities.xml
       rel                                           title templated
1     self                                 Mosaic metadata        NA
2 tilejson                   TileJSON link (Template URL).      TRUE
3 tilejson   TileJSON link for `agb` layer (Template URL).      TRUE
4      map                 Map viewer link (Template URL).      TRUE
5      map Map viewer link for `agb` layer (Template URL).      TRUE
6     wmts                        WMTS link (Template URL)      TRUE

Render Tiles for the Whole Collection

The /collections endpoint in a titiler-pgstac deployment enables rendering tiles for an entire STAC collection as a mosaic, using just the collection ID. This is useful when collections have full spatial coverage per item (e.g., for aggregated or max-reduced layers).

The code below queries the /collections/{collection_id}/info endpoint to retrieve default rendering parameters for the "agb" layer, including the assets (e.g., "tif") and the rescaling range. These are then passed to the TileJSON URL with the correct tile matrix ID (WebMercatorQuad) to get back a fully configured tile URL. The parsed TileJSON response includes the actual tile URL template along with min/max zoom, bounds, and metadata for client-side map rendering.

[4]:

agb_raw <- collection_info$search$metadata$defaults$agb

agb_render_params <- list(
  assets = agb_raw$assets,
  rescale = paste(agb_raw$rescale[1, ], collapse = ","),
  colormap_name = agb_raw$colormap_name,
  bidx = 1
)

print(agb_render_params)


tilejson_url <- collection_info$links$href[
  collection_info$links$rel == "tilejson"
][[1]]
tilejson_url <- gsub("\\{tileMatrixSetId\\}", "WebMercatorQuad", tilejson_url)

$assets
[1] "tif"

$rescale
[1] "0,400"

$colormap_name
[1] "gist_earth_r"

$bidx
[1] 1

[5]:

response_tilejson <- request(tilejson_url) |>
  req_url_query(!!!agb_render_params) |>
  req_perform()

collection_tilejson <- fromJSON(resp_body_string(response_tilejson), simplifyVector = TRUE)
print(collection_tilejson)

$tilejson
[1] "2.2.0"

$name
[1] "Mosaic for 'icesat2-boreal' Collection"

$version
[1] "1.0.0"

$scheme
[1] "xyz"

$tiles
[1] "https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/tiles/WebMercatorQuad/{z}/{x}/{y}?assets=tif&rescale=0%2C400&colormap_name=gist_earth_r&bidx=1"

$minzoom
[1] 0

$maxzoom
[1] 24

$bounds
[1] -180.0   51.6  180.0   78.0

$center
[1]  0.0 64.8  0.0

[6]:

tile_url <- collection_tilejson$tiles[[1]]
min_zoom <- collection_tilejson$minzoom
max_zoom <- collection_tilejson$maxzoom

Leaflet Map with AGB Tiles and Bounding Box

The bounding box from collection_info$search$metadata$bounds is converted to an sf polygon. A Leaflet map is rendered with OpenStreetMap as base and AGB raster tiles (from TileJSON) as overlay. The tile layer uses minzoom/maxzoom parameters and is clipped to valid longitudes. A centroid-based setView() centers the map, and interactive layers are toggleable via controls.

[7]:

bbox_vals <- collection_info$search$metadata$bounds
names(bbox_vals) <- c("xmin", "ymin", "xmax", "ymax")

bbox_vals["xmin"] <- max(bbox_vals["xmin"], -179.9)
bbox_vals["xmax"] <- min(bbox_vals["xmax"], 179.9)

bbox <- st_bbox(bbox_vals, crs = st_crs(4326))
bbox_sf <- st_as_sfc(bbox) |> st_sf(geometry = _)

[8]:

center_lng <- collection_tilejson$center[1]
center_lat <- collection_tilejson$center[2]

leaflet() %>%
  addProviderTiles("OpenStreetMap", group = "Base Map") %>%
  addTiles(
    urlTemplate = tile_url,
    options = tileOptions(minZoom = min_zoom, maxZoom = max_zoom, opacity = 1),
    group = paste(collection_id, "- agb")
  ) %>%
  addPolygons(
    data = bbox_sf,
    color = "orange",
    weight = 2,
    opacity = 1,
    fillOpacity = 0,
    group = "Collection bounds"
  ) %>%
  addLayersControl(
    overlayGroups = c(paste(collection_id, "- agb"), "Collection bounds"),
    baseGroups = "Base Map",
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  setView(
    lng = center_lng,
    lat = center_lat,
    zoom = 5
  )

Render tiles for a STAC search

This example sends a POST request to the /searches/register endpoint on Titiler-PgSTAC. In the example below, we register a STAC search for the "icesat2-boreal" collection limited to the bounding box of Alaska. The request body includes collections and bbox as JSON.

The response contains a unique search ID and links to downstream tile rendering endpoints.

[9]:

# Alaska bounding box and collection ID
ak_bbox <- c(-180, 50, -128, 73)
collection_id <- "icesat2-boreal"


search_payload <- list(
  collections = list(collection_id),
  bbox = ak_bbox
)


response <- request("https://titiler-pgstac.maap-project.org/searches/register") |>
  req_body_json(search_payload) |>
  req_perform()


ak_search_request <- fromJSON(resp_body_string(response), flatten = TRUE)
print(ak_search_request)

$id
[1] "1bf89df0c73f2789c69e3c8ae3e5549f"

$links
                                                                                                                      href
1                                   https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/info
2        https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/{tileMatrixSetId}/tilejson.json
3             https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/{tileMatrixSetId}/map.html
4 https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/{tileMatrixSetId}/WMTSCapabilities.xml
       rel                              title templated
1 metadata                    Mosaic metadata        NA
2 tilejson   Link for TileJSON (Template URL)      TRUE
3      map Link for Map viewer (Template URL)      TRUE
4     wmts       Link for WMTS (Template URL)      TRUE

Retrieve TileJSON Using Search ID

This step constructs a tilejson URL using the search_id from the STAC query registration. To ensure proper styling, it manually includes the agb_render_params (e.g., rescale). A GET request is then sent to the constructed endpoint, and the JSON response is parsed. The result provides the visualization-ready tile URL and metadata for rendering.

[10]:

search_id <- ak_search_request[["id"]]
search_tilejson_url <- paste0(
  "https://titiler-pgstac.maap-project.org/searches/", search_id, "/WebMercatorQuad/tilejson.json"
)

agb_render_params$rescale <- paste(agb_render_params$rescale, collapse = ",")

response <- request(search_tilejson_url) |>
  req_url_query(!!!agb_render_params) |>
  req_perform()

ak_tilejson <- fromJSON(resp_body_string(response), simplifyVector = TRUE)
print(ak_tilejson)

$tilejson
[1] "2.2.0"

$name
[1] "1bf89df0c73f2789c69e3c8ae3e5549f"

$version
[1] "1.0.0"

$scheme
[1] "xyz"

$tiles
[1] "https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/tiles/WebMercatorQuad/{z}/{x}/{y}?assets=tif&rescale=0%2C400&colormap_name=gist_earth_r&bidx=1"

$minzoom
[1] 0

$maxzoom
[1] 24

$bounds
[1] -180   50 -128   73

$center
[1] -154.0   61.5    0.0

Render Leaflet Map from STAC Search Results

This block creates a Leaflet map in R that visualizes GEDI data tiles filtered by a STAC search. It extracts the tile rendering URL from the tilejson link and replaces {tileMatrixSetId} with "WebMercatorQuad". Rendering parameters such as rescale are added manually to the request. The response is parsed into a tilejson object and used to overlay GEDI tiles and STAC search bounds on the map.

[11]:

center_lng <- ak_tilejson$center[1]
center_lat <- ak_tilejson$center[2]

leaflet() %>%
  addProviderTiles("OpenStreetMap", group = "Base Map") %>%
  addTiles(
    urlTemplate = ak_tilejson$tiles[[1]],
    options = tileOptions(
      minZoom = ak_tilejson$minzoom,
      maxZoom = ak_tilejson$maxzoom,
      opacity = 1
    ),
    group = "Search Tiles"
  ) %>%
  addPolygons(
    data = bbox_sf,
    color = "orange",
    weight = 2,
    opacity = 1,
    fillOpacity = 0,
    group = "STAC search bounds"
  ) %>%
  addLayersControl(
    overlayGroups = c("Search Tiles", "STAC search bounds"),
    baseGroups = "Base Map",
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  setView(lng = center_lng, lat = center_lat, zoom = 6)