Visualizing rasters with MosaicJSON using stac_ipyleaflet

%pip install cogeo_mosaic --quiet

Run This Notebook

To access and run this tutorial within MAAP’s Algorithm Development Environment (ADE), please refer to the “Getting started with the MAAP” section of our documentation.

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.

# Set variables and helper functions
import os

titiler_url = "https://titiler.maap-project.org"
min_zoom = 0
max_zoom = 0
band_min = 0
band_max = 4000
user = os.getenv('CHE_WORKSPACE_NAMESPACE')
bucket = "maap-ops-workspace"
color_map = "gist_earth_r"

def checkFilePath(file_path):
    result = s3.list_objects(Bucket=bucket, Prefix=file_path)
    exists = True if 'Contents' in result else False
    if exists:
        print('PATH EXISTS')
        return [i for i in result['Contents'] if i["Key"].endswith('.tif')]
    return exists

Find data location

import boto3

path = "/local/path/to/files/" # "Path to directory with rasters:", Right-click on the directory and select option to Copy Path

s3 = boto3.client('s3')
file_name = path.split('/', 1)[-1]
if 'shared-buckets' in path:
    file_path = f'shared/{file_name}'
if 'my-private-bucket' in path:
    file_path = f'{user}/{file_name}'
if 'my-public-bucket' in path:
    file_path = f'shared/{user}/{file_name}'


# If using your files
#if file_path:
#    paths = checkFilePath(file_path)
#    files = [ i['Key'] for i in paths ]

# Our sample files
files = ['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']

print(files)

['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']

Get raster info (bounds, zoom, data type) from first file in array of files

import httpx

url = files[0]

r = httpx.get(
    f"{titiler_url}/cog/info",
    params = {
        "url": url,
    }
).json()

if r.get('count') > 0:
    bounds = r.get("bounds")
    min_zoom = r.get("minzoom")
    max_zoom = r.get("maxzoom")
    zoom = min_zoom + 1 if min_zoom == 0 else min_zoom
    bands = r.get("band_metadata")

    print("Bounds:", bounds)
    print("Zoom:", zoom, "min =", min_zoom, "max =", max_zoom)
    print("Data type:", r.get("dtype"))
    print("Bands:", bands)

Bounds: [-117.10749852280769, 50.78795362739066, -116.50936927974429, 51.16389512140189]
Zoom: 9 min = 9 max = 12
Data type: float32
Bands: [['b1', {}], ['b2', {}], ['b3', {}], ['b4', {}], ['b5', {}], ['b6', {}], ['b7', {}], ['b8', {}], ['b9', {}], ['b10', {}], ['b11', {}], ['b12', {}], ['b13', {}], ['b14', {}], ['b15', {}], ['b16', {}], ['b17', {}], ['b18', {}]]

Create mosaic from file urls

from cogeo_mosaic.mosaic import MosaicJSON

if min_zoom == 0 and max_zoom == 0:
    print("Warning: missing zoom attributes!")

mosaicdata = MosaicJSON.from_urls(files, minzoom=min_zoom, maxzoom=max_zoom)
print("Mosaic created!")

Mosaic created!

Upload the mosaicjson to the TiTiler and filter for tilejson results

mosaic_links = httpx.post(
    url=f"{titiler_url}/mosaics",
    headers={
        "Content-Type": "application/vnd.titiler.mosaicjson+json",
    },
    json=mosaicdata.model_dump(exclude_none=True),
).json()
# print(mosaic_links)

tilejson_object = list(
    filter(lambda x: x.get("rel") == "tilejson", dict(mosaic_links)["links"])
)
tilejson_url = tilejson_object[0]["href"]
tilejson_url

'https://titiler.maap-project.org/mosaics/af2ee61a-b8ee-4ee9-afcc-98b4f4b13ce3/tilejson.json'

Calculate raster center for map placement

r = httpx.get(tilejson_url).json()
center_data = r.get('center')
center = (center_data[1], center_data[0])
zoom = center_data[2]

from shapely.geometry import box

polygon = box(*bounds)
center = (polygon.centroid.y, polygon.centroid.x)
print("Center:", center)

Center: (50.97592437439628, -116.80843390127599)

Create the TileLayer

from ipyleaflet import TileLayer

params = {
    "return_mask": "true",
    "rescale": f"{band_min}, {band_max}",
    "bidx": "1",
    "colormap_name": "viridis"
}

r = httpx.get(tilejson_url, params=params).json()

layer_url = r['tiles'][0]
custom_layer = TileLayer(url=layer_url, show_loading=True, transparent=True)

Add the mosaic tile layer to the map

import stac_ipyleaflet

m = stac_ipyleaflet.StacIpyleaflet(zoom=zoom, center=center)
m.add_layer(custom_layer)
m