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.

[1]:
import glob
import os
import urllib

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.

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

# 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/maap-users")

tiles = [local_to_s3(file) for file in files]
print(tiles)
['s3://maap-ops-workspace/maap-users/alexdevseed/landsat8/viz/Landsat8_30542_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/maap-users/alexdevseed/landsat8/viz/Landsat8_30543_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/maap-users/alexdevseed/landsat8/viz/Landsat8_30822_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/maap-users/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(files, minzoom=9, maxzoom=16)

[4]:
# Now save the MosaicJSON to a file on S3 so the tile server can use it.
with MosaicBackend("L8_tile_test.json", mosaic_def=mosaicdata) as mosaic:
    mosaic.write(overwrite=True)

# The S3 path to the file is going to be
# s3://maap-ops-workspace/maap-users/alexdevseed/L8_tile_test.json

Make a Map

[5]:
import ipycmc
[6]:
w = ipycmc.MapCMC()
w

ADE Mosaic

[7]:
# 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 = "https://titiler.maap-project.org/mosaicjson/WMTSCapabilities.xml"
params = {
    "tile_format":"png",
    "tile_scale":"1",
    "pixel_selection":"first",
    "TileMatrixSetId":"WebMercatorQuad",
    "url":"s3://maap-ops-workspace/maap-users/alexdevseed/landsat8/viz/L8_tile_test.json",
    "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

[7]:
'https://h9su0upami.execute-api.us-east-1.amazonaws.com//mosaicjson/WMTSCapabilities.xml?tile_format=png&tile_scale=1&pixel_selection=first&TileMatrixSetId=WebMercatorQuad&url=s3%3A%2F%2Fmaap-ops-workspace%2Fmaap-users%2Falexdevseed%2Flandsat8%2Fviz%2FL8_tile_test.json&bidx=6&resampling_method=nearest&rescale=0%2C1&return_mask=true&colormap_name=viridis'
[8]:
# This adds a new layer to the map above, call Cloud Optimized GeoTIFF
w.load_layer_config(wmts_call, "wmts/xml")