{ "cells": [ { "cell_type": "markdown", "id": "5a5f9b7a-2bd9-40e9-97be-e885773ffa09", "metadata": {}, "source": [ "# Visualizing MAAP STAC Dataset with titiler-pgstac in R" ] }, { "cell_type": "markdown", "id": "aae2acc8-38fc-45da-9391-7ab12d8a70a0", "metadata": {}, "source": [ "Authors: Harshini Girish (UAH), Sheyenne Kirkland (UAH), Alex Mandel (Development Seed), Henry Rodman (Development Seed)\n", "\n", "Date: July 16, 2025\n", "\n", "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](https://titiler-pgstac.maap-project.org/).\n" ] }, { "cell_type": "markdown", "id": "d96bf134-2f78-4468-b416-f98cf63a6b2b", "metadata": {}, "source": [ "## Run This Notebook\n", "\n", "To access and run this tutorial within MAAP's Algorithm Development Environment (ADE), please refer to the [\"Getting started with the MAAP\"](https://docs.maap-project.org/en/latest/getting_started/getting_started.html) section of our documentation.\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "93c75a2f-ee1b-4c06-aede-469261744de9", "metadata": { "jp-MarkdownHeadingCollapsed": true }, "source": [ "## Additional Resources\n", "\n", "- [Visualizing STAC Items](https://docs.maap-project.org/en/latest/technical_tutorials/visualizing.html): Official MAAP tutorial demonstrating how to visualize spatial footprints from STAC collections using different functions inside the ADE.\n", "\n", "- [Titiler-PgSTAC API Reference](https://titiler-pgstac.maap-project.org/api.html): Full API reference for interacting with the Titiler-PgSTAC service, which enables dynamic tiling and rendering of geospatial assets from MAAP's STAC catalog.\n", "\n", "- [Visualizing with Titiler-PgSTAC](https://docs.maap-project.org/en/latest/technical_tutorials/visualization/visualizing_titiler-pgstac.html): MAAP tutorial on rendering and visualizing COGs using the Titiler-PgSTAC tiling service, including usage examples and parameter explanations.\n" ] }, { "cell_type": "markdown", "id": "50b05067-97e4-4578-b79d-7bc208b8da9c", "metadata": {}, "source": [ "## Install/Import Libraries\n", "Let’s install and load the packages necessary for this tutorial." ] }, { "cell_type": "code", "execution_count": 1, "id": "e6ec781f-fb50-4088-a44d-77f5b7b7ec6d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE\n", "\n" ] } ], "source": [ "library(httr2)\n", "library(jsonlite)\n", "library(leaflet)\n", "library(rstac)\n", "library(sf)" ] }, { "cell_type": "markdown", "id": "fafde98e-83e6-4556-a9b6-dc8b63c17aab", "metadata": {}, "source": [ "## Fetch Collection Info from Endpoint\n", "\n", "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. \n", "This information is used to configure how the collection is rendered as a map layer.\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "b1afbaf1-4727-4690-a521-da9effd0066b", "metadata": {}, "outputs": [], "source": [ "titiler_endpoint <- \"https://titiler-pgstac.maap-project.org\"\n", "collection_id <- \"icesat2-boreal\"\n", "info_url <- paste0(titiler_endpoint, \"/collections/\", collection_id, \"/info\")\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "c2fd72b2-b4df-49d6-8e92-3df2037b013d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "$search\n", "$search$hash\n", "[1] \"e1458a47fc19d7e34c06f3d8469fd4f6\"\n", "\n", "$search$search\n", "$search$search$collections\n", "[1] \"icesat2-boreal\"\n", "\n", "$search$search$`filter-lang`\n", "[1] \"cql2-json\"\n", "\n", "\n", "$search$lastused\n", "[1] \"2025-12-18T20:40:25.515213Z\"\n", "\n", "$search$usecount\n", "[1] 175\n", "\n", "$search$metadata\n", "$search$metadata$type\n", "[1] \"mosaic\"\n", "\n", "$search$metadata$bounds\n", "[1] -180.0 51.6 180.0 78.0\n", "\n", "$search$metadata$name\n", "[1] \"Mosaic for 'icesat2-boreal' Collection\"\n", "\n", "$search$metadata$assets\n", "[1] \"csv\" \"tif\"\n", "\n", "$search$metadata$defaults\n", "$search$metadata$defaults$agb\n", "$search$metadata$defaults$agb$bidx\n", "[1] 1\n", "\n", "$search$metadata$defaults$agb$assets\n", "[1] \"tif\"\n", "\n", "$search$metadata$defaults$agb$nodata\n", "[1] \"nan\"\n", "\n", "$search$metadata$defaults$agb$rescale\n", " [,1] [,2]\n", "[1,] 0 400\n", "\n", "$search$metadata$defaults$agb$colormap_name\n", "[1] \"gist_earth_r\"\n", "\n", "\n", "\n", "\n", "\n", "$links\n", " href\n", "1 https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/info\n", "2 https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/tilejson.json\n", "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\n", "4 https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/map.html\n", "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\n", "6 https://titiler-pgstac.maap-project.org/collections/icesat2-boreal/{tileMatrixSetId}/WMTSCapabilities.xml\n", " rel title templated\n", "1 self Mosaic metadata NA\n", "2 tilejson TileJSON link (Template URL). TRUE\n", "3 tilejson TileJSON link for `agb` layer (Template URL). TRUE\n", "4 map Map viewer link (Template URL). TRUE\n", "5 map Map viewer link for `agb` layer (Template URL). TRUE\n", "6 wmts WMTS link (Template URL) TRUE\n", "\n" ] } ], "source": [ "response <- request(info_url) |> req_timeout(10) |> req_perform()\n", "collection_info <- fromJSON(resp_body_string(response), simplifyVector = TRUE)\n", "print(collection_info)" ] }, { "cell_type": "markdown", "id": "bb59a6ea-4140-4a77-8a47-cbfa4afd7da4", "metadata": {}, "source": [ "## Render Tiles for the Whole Collection" ] }, { "cell_type": "markdown", "id": "f2e09482-ceb9-4b22-a5b5-4c98787b67e5", "metadata": {}, "source": [ "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).\n", "\n", "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.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "885863e6-96a7-4996-9b35-cc1622ba4468", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "$assets\n", "[1] \"tif\"\n", "\n", "$rescale\n", "[1] \"0,400\"\n", "\n", "$colormap_name\n", "[1] \"gist_earth_r\"\n", "\n", "$bidx\n", "[1] 1\n", "\n" ] } ], "source": [ "\n", "agb_raw <- collection_info$search$metadata$defaults$agb\n", "\n", "agb_render_params <- list(\n", " assets = agb_raw$assets,\n", " rescale = paste(agb_raw$rescale[1, ], collapse = \",\"),\n", " colormap_name = agb_raw$colormap_name,\n", " bidx = 1\n", ")\n", "\n", "print(agb_render_params)\n", "\n", "\n", "tilejson_url <- collection_info$links$href[\n", " collection_info$links$rel == \"tilejson\"\n", "][[1]]\n", "tilejson_url <- gsub(\"\\\\{tileMatrixSetId\\\\}\", \"WebMercatorQuad\", tilejson_url)\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "e6c0f9c7-a470-4050-bbca-143e4c7cc594", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "$tilejson\n", "[1] \"2.2.0\"\n", "\n", "$name\n", "[1] \"Mosaic for 'icesat2-boreal' Collection\"\n", "\n", "$version\n", "[1] \"1.0.0\"\n", "\n", "$scheme\n", "[1] \"xyz\"\n", "\n", "$tiles\n", "[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\"\n", "\n", "$minzoom\n", "[1] 0\n", "\n", "$maxzoom\n", "[1] 24\n", "\n", "$bounds\n", "[1] -180.0 51.6 180.0 78.0\n", "\n", "$center\n", "[1] 0.0 64.8 0.0\n", "\n" ] } ], "source": [ "response_tilejson <- request(tilejson_url) |>\n", " req_url_query(!!!agb_render_params) |>\n", " req_perform()\n", "\n", "collection_tilejson <- fromJSON(resp_body_string(response_tilejson), simplifyVector = TRUE)\n", "print(collection_tilejson)" ] }, { "cell_type": "code", "execution_count": 6, "id": "a650c338-88fe-4c3f-9256-270228b76a91", "metadata": {}, "outputs": [], "source": [ "tile_url <- collection_tilejson$tiles[[1]]\n", "min_zoom <- collection_tilejson$minzoom\n", "max_zoom <- collection_tilejson$maxzoom\n" ] }, { "cell_type": "markdown", "id": "6ded02fa-116b-4d66-8b85-4525a93f1454", "metadata": {}, "source": [ "## Leaflet Map with AGB Tiles and Bounding Box\n", "\n", "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.\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "08c9fc0c-0067-4eb2-bc20-66a0ab131acd", "metadata": {}, "outputs": [], "source": [ "bbox_vals <- collection_info$search$metadata$bounds\n", "names(bbox_vals) <- c(\"xmin\", \"ymin\", \"xmax\", \"ymax\")\n", "\n", "bbox_vals[\"xmin\"] <- max(bbox_vals[\"xmin\"], -179.9)\n", "bbox_vals[\"xmax\"] <- min(bbox_vals[\"xmax\"], 179.9)\n", "\n", "bbox <- st_bbox(bbox_vals, crs = st_crs(4326))\n", "bbox_sf <- st_as_sfc(bbox) |> st_sf(geometry = _)\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "8910470a-fe3f-4941-bd1a-6275663595cf", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\t\n", "\t\t\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\t\n", "\t\n", "\t\t
\n", "\n", "\t\n", "\n" ], "text/plain": [ "HTML widgets cannot be represented in plain text (need html)" ] }, "metadata": { "text/html": { "isolated": true } }, "output_type": "display_data" } ], "source": [ "center_lng <- collection_tilejson$center[1]\n", "center_lat <- collection_tilejson$center[2]\n", "\n", "leaflet() %>%\n", " addProviderTiles(\"OpenStreetMap\", group = \"Base Map\") %>%\n", " addTiles(\n", " urlTemplate = tile_url,\n", " options = tileOptions(minZoom = min_zoom, maxZoom = max_zoom, opacity = 1),\n", " group = paste(collection_id, \"- agb\")\n", " ) %>%\n", " addPolygons(\n", " data = bbox_sf,\n", " color = \"orange\",\n", " weight = 2,\n", " opacity = 1,\n", " fillOpacity = 0,\n", " group = \"Collection bounds\"\n", " ) %>%\n", " addLayersControl(\n", " overlayGroups = c(paste(collection_id, \"- agb\"), \"Collection bounds\"),\n", " baseGroups = \"Base Map\",\n", " options = layersControlOptions(collapsed = FALSE)\n", " ) %>%\n", " setView(\n", " lng = center_lng,\n", " lat = center_lat,\n", " zoom = 5\n", " )\n" ] }, { "cell_type": "markdown", "id": "8ed1fe28-d2e0-4562-a03c-0caaa05374eb", "metadata": {}, "source": [ "## Render tiles for a STAC search" ] }, { "cell_type": "markdown", "id": "d1664b83-8bd1-4f56-9886-bae7bad59fb0", "metadata": {}, "source": [ "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. \n", "The response contains a unique search ID and links to downstream tile rendering endpoints.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "9d8cbd30-0131-4047-8c71-163cc920b803", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "$id\n", "[1] \"1bf89df0c73f2789c69e3c8ae3e5549f\"\n", "\n", "$links\n", " href\n", "1 https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/info\n", "2 https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/{tileMatrixSetId}/tilejson.json\n", "3 https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/{tileMatrixSetId}/map.html\n", "4 https://titiler-pgstac.maap-project.org/searches/1bf89df0c73f2789c69e3c8ae3e5549f/{tileMatrixSetId}/WMTSCapabilities.xml\n", " rel title templated\n", "1 metadata Mosaic metadata NA\n", "2 tilejson Link for TileJSON (Template URL) TRUE\n", "3 map Link for Map viewer (Template URL) TRUE\n", "4 wmts Link for WMTS (Template URL) TRUE\n", "\n" ] } ], "source": [ "# Alaska bounding box and collection ID\n", "ak_bbox <- c(-180, 50, -128, 73)\n", "collection_id <- \"icesat2-boreal\"\n", "\n", "\n", "search_payload <- list(\n", " collections = list(collection_id),\n", " bbox = ak_bbox\n", ")\n", "\n", "\n", "response <- request(\"https://titiler-pgstac.maap-project.org/searches/register\") |>\n", " req_body_json(search_payload) |>\n", " req_perform()\n", "\n", "\n", "ak_search_request <- fromJSON(resp_body_string(response), flatten = TRUE)\n", "print(ak_search_request)\n" ] }, { "cell_type": "markdown", "id": "a6a49388-d7da-42be-9a57-bd00cd4e5d17", "metadata": {}, "source": [ "## Retrieve TileJSON Using Search ID\n", "\n", "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.\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "604ab20f-a30f-4da2-a8f8-087613d54844", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "$tilejson\n", "[1] \"2.2.0\"\n", "\n", "$name\n", "[1] \"1bf89df0c73f2789c69e3c8ae3e5549f\"\n", "\n", "$version\n", "[1] \"1.0.0\"\n", "\n", "$scheme\n", "[1] \"xyz\"\n", "\n", "$tiles\n", "[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\"\n", "\n", "$minzoom\n", "[1] 0\n", "\n", "$maxzoom\n", "[1] 24\n", "\n", "$bounds\n", "[1] -180 50 -128 73\n", "\n", "$center\n", "[1] -154.0 61.5 0.0\n", "\n" ] } ], "source": [ "search_id <- ak_search_request[[\"id\"]]\n", "search_tilejson_url <- paste0(\n", " \"https://titiler-pgstac.maap-project.org/searches/\", search_id, \"/WebMercatorQuad/tilejson.json\"\n", ")\n", "\n", "agb_render_params$rescale <- paste(agb_render_params$rescale, collapse = \",\")\n", "\n", "response <- request(search_tilejson_url) |>\n", " req_url_query(!!!agb_render_params) |>\n", " req_perform()\n", "\n", "ak_tilejson <- fromJSON(resp_body_string(response), simplifyVector = TRUE)\n", "print(ak_tilejson)\n" ] }, { "cell_type": "markdown", "id": "0ff735e3-6ceb-417b-ac1c-099999e8d3ea", "metadata": {}, "source": [ "## Render Leaflet Map from STAC Search Results\n", "\n", "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.\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "c8b75c0d-1a80-4c70-9558-3877fed362cd", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\t\n", "\t\t\n", "\t\t\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\t\n", "\t\n", "\t\t
\n", "\n", "\t\n", "\n" ], "text/plain": [ "HTML widgets cannot be represented in plain text (need html)" ] }, "metadata": { "text/html": { "isolated": true } }, "output_type": "display_data" } ], "source": [ "center_lng <- ak_tilejson$center[1]\n", "center_lat <- ak_tilejson$center[2]\n", "\n", "leaflet() %>%\n", " addProviderTiles(\"OpenStreetMap\", group = \"Base Map\") %>%\n", " addTiles(\n", " urlTemplate = ak_tilejson$tiles[[1]], \n", " options = tileOptions(\n", " minZoom = ak_tilejson$minzoom,\n", " maxZoom = ak_tilejson$maxzoom,\n", " opacity = 1\n", " ),\n", " group = \"Search Tiles\"\n", " ) %>%\n", " addPolygons(\n", " data = bbox_sf,\n", " color = \"orange\",\n", " weight = 2,\n", " opacity = 1,\n", " fillOpacity = 0,\n", " group = \"STAC search bounds\"\n", " ) %>%\n", " addLayersControl(\n", " overlayGroups = c(\"Search Tiles\", \"STAC search bounds\"),\n", " baseGroups = \"Base Map\",\n", " options = layersControlOptions(collapsed = FALSE)\n", " ) %>%\n", " setView(lng = center_lng, lat = center_lat, zoom = 6)\n" ] } ], "metadata": { "kernelspec": { "display_name": "R 4.4.3", "language": "R", "name": "ir4" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "4.4.3" } }, "nbformat": 4, "nbformat_minor": 5 }