Algorithm Development Environment (ADE) Visualization Example

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

Visualizing COGs generated in the MAAP Algorithm Development Environment (ADE)

This notebook visualizes COGs generated in the ADE using the pythonCommon Mapping Client, an open source project of NASA and JPL.

Any Cloud Optimized GeoTIFF (COG), or group of COGs, in an ADE workspace can be visualized in a dynamic map by using a tiling service hosted in MAAP.

Steps: 1. Make a list of TIFFs in your workspace to use as a single layer 2. Generate a MosaicJSON file from this list of files (or a GeoJSON index) 3. Combine the MosaicJSON with other tiler visualization parameters to register a layer with your visualization tool.

[11]:
import glob
import ipycmc
import os
from pprint import pprint
import requests
import urllib

#!pip install cogeo_mosaic
from cogeo_mosaic.mosaic import MosaicJSON
from cogeo_mosaic.backends import MosaicBackend

Build a list of files

You can either make a list of file paths, or create a GeoJSON layer with a column containing the file paths. The paths need to be S3 paths currently.

[12]:
# Local Path to your COGs
dps_output = "/projects/shared-buckets/alexdevseed/landsat8/viz/"

# titiler endpoint
titiler_endpoint = f"https://titiler.maap-project.org"

# Search for files to include, use recursive if nested folders (common in DPS output)
files = glob.glob(os.path.join(dps_output, "Landsat*_dps.tif"), recursive=False)

def local_to_s3(url):
    ''' A Function to convert local paths to s3 urls'''
    return url.replace("/projects/shared-buckets", "s3://maap-ops-workspace/shared")

tiles = [local_to_s3(file) for file in files]
print(tiles)
['s3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30542_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30543_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30822_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30823_comp_cog_2015-2020_dps.tif']

How to find the S3 path

You might be wondering how to find the S3 path to begin with. Right now the easiest way is to right click on a file in the file explorer on the left panel, and Get Presigned S3 Url.

It will look something like https://maap-ops-workspace.s3.amazonaws.com/shared/alexdevseed/landsat8/viz/Copernicus_30438_covars_cog_topo_stack.tif?AWSAccessKeyId...

The first part of the url is the bucket name: maap-ops-workspace. After the next /, it then matches to the local path.

Future versions of MAAP will include functions to do this part for you…

Make a mosaic

[3]:
# Now take the list of S3 paths and generate a mosaic
# TODO: if you have a lot of files (more than 100), creating a GeoJSON index and using from_features will be more efficient.

mosaicdata = MosaicJSON.from_urls(tiles, minzoom=9, maxzoom=16)
[13]:
r = requests.post(
    url=f"{titiler_endpoint}/mosaics",
    headers={
        "Content-Type": "application/vnd.titiler.mosaicjson+json",
    },
    json=mosaicdata.dict(exclude_none=True)).json()

mosaicid = r['id']

Make a Map

[14]:
w = ipycmc.MapCMC()
w

ADE Mosaic

[18]:
# Build a WMTS call
"""
All of this is subject to change in a future version
The important parameters for users:
  url : the S3 path to the MosaicJSON file,
  bidx (band number),
  rescale (required if using non Byte data type),
  colormap_name or colormap

Other parameters are possible, see https://titiler.maap-project.org/docs#/MosaicJSON/wmts_mosaicjson_WMTSCapabilities_xml_get
"""

wmts_url = f"https://titiler.maap-project.org/mosaics/{mosaicid}/WMTSCapabilities.xml"
params = {
    "tile_format":"png",
    "tile_scale":"1",
    "pixel_selection":"first",
    "TileMatrixSetId":"WebMercatorQuad",
    "bidx":"6", # Select which band to use
    "resampling_method":"nearest",
    "rescale":"0,1", # Values in data are from 0 to 1
    "return_mask":"true",
    "colormap_name":"viridis" # Any colormap from matplotlib will work
}

wmts_call = "?".join([wmts_url, urllib.parse.urlencode(params)])

# Note Jupyter bug add amp; incorrectly when printing the url
wmts_call

[18]:
'https://titiler.maap-project.org/mosaics/d96e662d-fa69-4de1-a36c-4cf144c3e8e9/WMTSCapabilities.xml?tile_format=png&tile_scale=1&pixel_selection=first&TileMatrixSetId=WebMercatorQuad&bidx=6&resampling_method=nearest&rescale=0%2C1&return_mask=true&colormap_name=viridis'
[16]:
# This adds a new layer to the map above, call Cloud Optimized GeoTIFF
w.load_layer_config(wmts_call, "wmts/xml")