Visualizing SRTM dataset from MAAP CMR STAC using MosaicJSON

In this notebook, we discover SRTM COGs from `MAAP’s CMR STAC API Proxy <http://cmr-stac-listener-lb-uat-807915637.us-east-1.elb.amazonaws.com/stac/NASA_MAAP/collections/SRTMGL1_COD.v001>`__ and generate a mosaic.

MosaicJSON

A common challenge in visualizing spatial data is data is stored across many files representing small spatial extents.

MosaicJSON is a specification created by DevelopmentSeed which aims to be an open standard for representing metadata about a spatial mosaic of many COG files.

MosaicJSON can be seen as a cloud friendly virtual raster (see GDAL’s VRT) enabling spatial and temporal indexing for a list of Cloud-Optimized GeoTIFF.

Ref: https://github.com/developmentseed/mosaicjson-spec

Data

The NASA Shuttle Radar Topographic Mission (SRTM) has provided digital elevation data (DEMs) for over 80% of the globe. This data is currently distributed free of charge by USGS and is available for download from the National Map Seamless Data Distribution System, or the USGS FTP site.

At MAAP, we’ve converted this elevation data into Cloud-Optimized GeoTIFFs (COGs) so they can be efficiently queried and visualized. These COGs are available in the MAAP CMR.

Titiler Endpoint

By default, TiTiler has a complete mosaicjson endpoint. For this demo we are going to use an instance of TiTiler deployed by MAAP at https://h9su0upami.execute-api.us-east-1.amazonaws.com/

Requirements

To be able to run this notebook you’ll need the following requirements: - rasterio - folium - requests - tqdm - rio-tiler (2.0b8) (Optional) - cogeo-mosaic (Optional)

pip install rasterio folium requests tqdm requests

pip install rio-tiler cogeo-mosaic --pre

Get the Data

[1]:
import os
import json
import rasterio
import requests
from pprint import pprint
from concurrent import futures
import urllib

from rio_tiler.io import COGReader
from rasterio.features import bounds as featureBounds

from folium import Map, TileLayer, GeoJson

Fetch SRTM COG STAC items

[2]:
stac_endpoint = "https://cmr-stac.dit.maap-project.org/stac/NASA_MAAP/search"

stac_json = {
    "collections": ["SRTMGL1_COD.v001"],
    "bbox": "4,42,16,48",
    "limit": 120
}

headers = {
            "Content-Type": "application/json",
            "Accept": "application/geo+json, application/json",
        }

# Read Items

r_stac = requests.post(stac_endpoint, headers=headers, json=stac_json)
items = r_stac.json()

Map the data bounds

[3]:
geojson = {'type': 'FeatureCollection', 'features': items['features']}

bounds = featureBounds(geojson)
zoom_start = 5

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

geo_json = GeoJson(
    data=geojson,
    style_function=lambda x: {
        'opacity': 1, 'dashArray': '1', 'fillOpacity': 0, 'weight': 1
    },
)
geo_json.add_to(m)
m
[3]:
Make this Notebook Trusted to load map: File -> Trust Notebook

2. Create Mosaic

We’re using the TiTiler deployed by MAAP

[4]:
titiler_endpoint = "https://h9su0upami.execute-api.us-east-1.amazonaws.com/"  # MAAP titiler endpoint
[5]:
%time
from rio_tiler.io import COGReader

first_cog = geojson['features'][0]['assets']['browse']['href']
with COGReader(first_cog) as cog:
    info = cog.info()
CPU times: user 1 µs, sys: 2 µs, total: 3 µs
Wall time: 6.91 µs
[6]:
from cogeo_mosaic.mosaic import MosaicJSON

# We are creating the mosaicJSON using the geojson
# SRTMGL1 CODs have a "browse" asset type, so we can create a mosaic of these type="image/tiff" assets
# accesor function provide access to those
mosaicdata = MosaicJSON.from_features(geojson.get('features'), minzoom=6, maxzoom=info.maxzoom, accessor=lambda feature : feature['assets']['browse']['href'])
[7]:
# Uploading the mosaicjson to the TiTiler
r = requests.post(
    url=f"{titiler_endpoint}/mosaicjson/mosaics",
    headers={
        "Content-Type": "application/vnd.titiler.mosaicjson+json",
    },
    json=mosaicdata.dict(exclude_none=True)).json()

pprint(r)
{'id': 'f53f5571-ee40-494e-a538-c49f56a9ffeb',
 'links': [{'href': 'https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/f53f5571-ee40-494e-a538-c49f56a9ffeb',
            'rel': 'self',
            'title': 'Self',
            'type': 'application/json'},
           {'href': 'https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/f53f5571-ee40-494e-a538-c49f56a9ffeb/mosaicjson',
            'rel': 'mosaicjson',
            'title': 'MosiacJSON',
            'type': 'application/json'},
           {'href': 'https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/f53f5571-ee40-494e-a538-c49f56a9ffeb/tilejson.json',
            'rel': 'tilejson',
            'title': 'TileJSON',
            'type': 'application/json'},
           {'href': 'https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/f53f5571-ee40-494e-a538-c49f56a9ffeb/tiles/{z}/{x}/{y}',
            'rel': 'tiles',
            'title': 'Tiles',
            'type': 'application/json'},
           {'href': 'https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/f53f5571-ee40-494e-a538-c49f56a9ffeb/WMTSCapabilities.xml',
            'rel': 'wmts',
            'title': 'WMTS',
            'type': 'application/json'}]}

The response from the post request gives endpoints for different services (eg. mosaicjson, tilejson, tiles, wmts, etc). We’re fetching the tilejson endpoint.

[8]:
tilejson_endpoint = list(filter(lambda x: x.get("rel") == "tilejson", dict(r)["links"]))

Display Tiles

NOTE: You have to zoom to “minzoom” to load the data.

SECOND NOTE: The important bit is the “tiles” endpoint returned from f"{titiler_endpoint}/mosaicjson/mosaics. This endpoint (e.g. https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/4199126b-9313-435a-b4b5-17802716b7b1/tiles/{z}/{x}/{y}) could be used in any map client which can tile x,y,z layers.

[9]:
r_te = requests.get(
    tilejson_endpoint[0].get('href')
).json()

pprint(r_te)

tiles = TileLayer(
    tiles=f"{r_te['tiles'][0]}rescale=0,4000",
    min_zoom=r_te["minzoom"],
    max_zoom=r_te["maxzoom"],
    opacity=1,
    attr="USGS"
)

tiles.add_to(m)
m
{'bounds': [2.9997222, 40.9997222, 17.0002778, 49.0002778],
 'center': [10.0, 45.0, 6],
 'maxzoom': 12,
 'minzoom': 6,
 'name': 'f53f5571-ee40-494e-a538-c49f56a9ffeb',
 'scheme': 'xyz',
 'tilejson': '2.2.0',
 'tiles': ['https://h9su0upami.execute-api.us-east-1.amazonaws.com/mosaicjson/mosaics/f53f5571-ee40-494e-a538-c49f56a9ffeb/tiles/{z}/{x}/{y}@1x?'],
 'version': '1.0.0'}
[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook