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 python `Common Mapping Client <https://github.com/nasa/common-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.

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.

# 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

# 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)

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

w = ipycmc.MapCMC()
w

# 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.ops.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

'https://titiler.ops.maap-project.org/mosaics/e613b570-7318-4d88-92c4-1c5373c6b513/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'

# This adds a new layer to the map above, call Cloud Optimized GeoTIFF
w.load_layer_config(wmts_call, "wmts/xml")