Working With MosaicJSON¶

MosaicJSON¶

MosaicJSON is a specification created by DevelopmentSeed which aims to be an open standard for representing metadata about a mosaic of Cloud-Optimized GeoTIFF (COG) files.

MosaicJSON can be seen as a Virtual raster (see GDAL's VRT) enabling spatial and temporal processing for a list of Cloud-Optimized GeoTIFF.

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

Data¶

For this demo, we are going to use CloudOptimized GeoTIFF from NOAA/Emergency Response Imagery: https://registry.opendata.aws/noaa-eri/

Requirement: AWS credentials

Endpoint¶

By default, TiTiler has mosaicjson endpoints.

Docs: https://titiler.xyz/api.html#/MosaicJSON

Requirements¶

To be able to run this notebook you'll need the following requirements:

• rasterio
• folium
• httpx
• tqdm
• rio-tiler
• cogeo-mosaic
• boto3
• geojson_pydantic

pip install rasterio folium tqdm httpx rio-tiler geojson_pydantic cogeo-mosaic

In [ ]:

# Uncomment this line if you need to install the dependencies
#!pip install rasterio folium tqdm httpx rio-tiler geojson_pydantic cogeo-mosaic

# Uncomment this line if you need to install the dependencies #!pip install rasterio folium tqdm httpx rio-tiler geojson_pydantic cogeo-mosaic

Get the Data¶

In [ ]:

import os
import json
import rasterio
import httpx
from boto3.session import Session as boto3_session

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

from folium import Map, TileLayer, GeoJson

import os import json import rasterio import httpx from boto3.session import Session as boto3_session from concurrent import futures from rio_tiler.io import COGReader from rasterio.features import bounds as featureBounds from folium import Map, TileLayer, GeoJson

1. Fetch and parse page¶

In [ ]:

# To Be able to run this notebook you'll need to have AWS credential available in the environment

# import os
# os.environ["AWS_ACCESS_KEY_ID"] = "YOUR AWS ACCESS ID HERE"
# os.environ["AWS_SECRET_ACCESS_KEY"] = "YOUR AWS ACCESS KEY HERE"

# To Be able to run this notebook you'll need to have AWS credential available in the environment # import os # os.environ["AWS_ACCESS_KEY_ID"] = "YOUR AWS ACCESS ID HERE" # os.environ["AWS_SECRET_ACCESS_KEY"] = "YOUR AWS ACCESS KEY HERE"

In [ ]:

session = boto3_session(region_name="us-west-2")
client = session.client("s3")

bucket = "noaa-eri-pds"  # https://registry.opendata.aws/omi-no2-nasa/


def list_objects(bucket, prefix):
    """AWS s3 list objects."""

    paginator = client.get_paginator("list_objects_v2")

    files = []
    for subset in paginator.paginate(Bucket=bucket, Prefix=prefix):
        files.extend(subset.get("Contents", []))

    return [r["Key"] for r in files]


files = list_objects(bucket, "2020_Nashville_Tornado/20200307a_RGB")
files = [f"s3://{bucket}/{f}" for f in files if f.endswith(".tif")]

print(f"Number of GeoTIFF: {len(files)}")

session = boto3_session(region_name="us-west-2") client = session.client("s3") bucket = "noaa-eri-pds" # https://registry.opendata.aws/omi-no2-nasa/ def list_objects(bucket, prefix): """AWS s3 list objects.""" paginator = client.get_paginator("list_objects_v2") files = [] for subset in paginator.paginate(Bucket=bucket, Prefix=prefix): files.extend(subset.get("Contents", [])) return [r["Key"] for r in files] files = list_objects(bucket, "2020_Nashville_Tornado/20200307a_RGB") files = [f"s3://{bucket}/{f}" for f in files if f.endswith(".tif")] print(f"Number of GeoTIFF: {len(files)}")

In [ ]:

print(files)

print(files)

2. Create Features and Viz (Optional)¶

Read each file geo metadata

In [ ]:

# We can derive the `bbox` from the filename
# s3://noaa-eri-pds/2020_Nashville_Tornado/20200307a_RGB/20200307aC0870130w361200n.tif
# -> 20200307aC0870130w361200n.tif
# -> 20200307aC "0870130w" "361200n" .tif
# -> 0870130w -> 87.025 (West)
# -> 361200n -> 36.2 (Top)
# We also know each files cover ~0.025x~0.025 degrees

import re
from geojson_pydantic.features import Feature
from geojson_pydantic.geometries import Polygon


def dms_to_degree(v: str) -> float:
    """convert degree minute second to decimal degrees.

    '0870130w' -> 87.025
    """
    deg = int(v[0:3])
    minutes = int(v[3:5])
    seconds = int(v[5:7])
    direction = v[-1].upper()
    return (float(deg) + float(minutes) / 60 + float(seconds) / (60 * 60)) * (
        -1 if direction in ["W", "S"] else 1
    )


def fname_to_feature(src_path: str) -> Feature:
    bname = os.path.basename(src_path)
    lon_dms = bname[10:18]
    lat_dms = bname[18:25]

    lon = dms_to_degree(lon_dms)
    lat = dms_to_degree("0" + lat_dms)

    return Feature(
        geometry=Polygon.from_bounds(lon, lat - 0.025, lon + 0.025, lat),
        properties={
            "path": src_path,
        },
    )


features = [fname_to_feature(f).dict(exclude_none=True) for f in files]

# OR We could use Rasterio/rio-tiler

# def worker(src_path: str) -> Feature:
#     try:
#         with COGReader(src_path) as cog:
#             wgs_bounds = cog.geographic_bounds
#     except:
#         return {}
#
#     return Feature(
#         geometry=Polygon.from_bounds(*wgs_bounds),
#         properties={
#             "path": src_path,
#         }
#     )
#
# with futures.ThreadPoolExecutor(max_workers=20) as executor:
#     features = [r.dict(exclude_none=True) for r in executor.map(worker, files) if r]

# We can derive the `bbox` from the filename # s3://noaa-eri-pds/2020_Nashville_Tornado/20200307a_RGB/20200307aC0870130w361200n.tif # -> 20200307aC0870130w361200n.tif # -> 20200307aC "0870130w" "361200n" .tif # -> 0870130w -> 87.025 (West) # -> 361200n -> 36.2 (Top) # We also know each files cover ~0.025x~0.025 degrees import re from geojson_pydantic.features import Feature from geojson_pydantic.geometries import Polygon def dms_to_degree(v: str) -> float: """convert degree minute second to decimal degrees. '0870130w' -> 87.025 """ deg = int(v[0:3]) minutes = int(v[3:5]) seconds = int(v[5:7]) direction = v[-1].upper() return (float(deg) + float(minutes) / 60 + float(seconds) / (60 * 60)) * ( -1 if direction in ["W", "S"] else 1 ) def fname_to_feature(src_path: str) -> Feature: bname = os.path.basename(src_path) lon_dms = bname[10:18] lat_dms = bname[18:25] lon = dms_to_degree(lon_dms) lat = dms_to_degree("0" + lat_dms) return Feature( geometry=Polygon.from_bounds(lon, lat - 0.025, lon + 0.025, lat), properties={ "path": src_path, }, ) features = [fname_to_feature(f).dict(exclude_none=True) for f in files] # OR We could use Rasterio/rio-tiler # def worker(src_path: str) -> Feature: # try: # with COGReader(src_path) as cog: # wgs_bounds = cog.geographic_bounds # except: # return {} # # return Feature( # geometry=Polygon.from_bounds(*wgs_bounds), # properties={ # "path": src_path, # } # ) # # with futures.ThreadPoolExecutor(max_workers=20) as executor: # features = [r.dict(exclude_none=True) for r in executor.map(worker, files) if r]

In [ ]:

geojson = {"type": "FeatureCollection", "features": features}

bounds = featureBounds(geojson)

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

geo_json = GeoJson(
    data=geojson,
    style_function=lambda x: {
        "opacity": 1,
        "dashArray": "1",
        "fillOpacity": 0,
        "weight": 1,
    },
)
geo_json.add_to(m)
m

geojson = {"type": "FeatureCollection", "features": features} bounds = featureBounds(geojson) m = Map( tiles="OpenStreetMap", location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom_start=6, ) geo_json = GeoJson( data=geojson, style_function=lambda x: { "opacity": 1, "dashArray": "1", "fillOpacity": 0, "weight": 1, }, ) geo_json.add_to(m) m

5. Create Mosaic¶

In [ ]:

from rio_tiler.io import COGReader
from cogeo_mosaic.mosaic import MosaicJSON
from cogeo_mosaic.backends import MosaicBackend

with COGReader(files[0]) as cog:
    info = cog.info()
    print(info.minzoom)
    print(info.maxzoom)

from rio_tiler.io import COGReader from cogeo_mosaic.mosaic import MosaicJSON from cogeo_mosaic.backends import MosaicBackend with COGReader(files[0]) as cog: info = cog.info() print(info.minzoom) print(info.maxzoom)

In [ ]:

# We are creating the mosaicJSON using the features we created earlier
# by default MosaicJSON.from_feature will look in feature.properties.path to get the path of the dataset
mosaicdata = MosaicJSON.from_features(
    features, minzoom=info.minzoom, maxzoom=info.maxzoom
)
with MosaicBackend("NOAA_Nashville_Tornado.json.gz", mosaic_def=mosaicdata) as mosaic:
    mosaic.write(overwrite=True)
    print(mosaic.info())

# We are creating the mosaicJSON using the features we created earlier # by default MosaicJSON.from_feature will look in feature.properties.path to get the path of the dataset mosaicdata = MosaicJSON.from_features( features, minzoom=info.minzoom, maxzoom=info.maxzoom ) with MosaicBackend("NOAA_Nashville_Tornado.json.gz", mosaic_def=mosaicdata) as mosaic: mosaic.write(overwrite=True) print(mosaic.info())

In [ ]:

titiler_endpoint = (
    "https://titiler.xyz"  # Developmentseed Demo endpoint. Please be kind.
)

r = httpx.get(
    f"{titiler_endpoint}/mosaicjson/WebMercatorQuad/tilejson.json",
    params={
        # For this demo we are use the same mosaic but stored on the web
        "url": "https://gist.githubusercontent.com/vincentsarago/c6ace3ccd29a82a4a5531693bbcd61fc/raw/e0d0174a64a9acd2fb820f2c65b1830aab80f52b/NOAA_Nashville_Tornado.json"
    },
).json()
print(r)

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

tiles = TileLayer(
    tiles=r["tiles"][0],
    min_zoom=r["minzoom"],
    max_zoom=r["maxzoom"],
    opacity=1,
    attr="NOAA",
)

geo_json = GeoJson(
    data=geojson,
    style_function=lambda x: {
        "opacity": 1,
        "dashArray": "1",
        "fillOpacity": 0,
        "weight": 1,
    },
)
tiles.add_to(m)
geo_json.add_to(m)
m

titiler_endpoint = ( "https://titiler.xyz" # Developmentseed Demo endpoint. Please be kind. ) r = httpx.get( f"{titiler_endpoint}/mosaicjson/WebMercatorQuad/tilejson.json", params={ # For this demo we are use the same mosaic but stored on the web "url": "https://gist.githubusercontent.com/vincentsarago/c6ace3ccd29a82a4a5531693bbcd61fc/raw/e0d0174a64a9acd2fb820f2c65b1830aab80f52b/NOAA_Nashville_Tornado.json" }, ).json() print(r) m = Map( location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom_start=13 ) tiles = TileLayer( tiles=r["tiles"][0], min_zoom=r["minzoom"], max_zoom=r["maxzoom"], opacity=1, attr="NOAA", ) geo_json = GeoJson( data=geojson, style_function=lambda x: { "opacity": 1, "dashArray": "1", "fillOpacity": 0, "weight": 1, }, ) tiles.add_to(m) geo_json.add_to(m) m

In [ ]: