stac_to_mosaicJSON.py

import mercantile
from shapely.geometry import box, shape
from supermercado import burntiles

minzoom = 7
maxzoom = 12
maximum_items_per_tile = 20

def stac_to_mosaicJSON(data):
    dataset = data["features"]

    tiles = burntiles.burn(dataset, minzoom)
    tiles = list(set(["{2}-{0}-{1}".format(*tile.tolist()) for tile in tiles]))
    bounds = burntiles.find_extrema(dataset)

    mosaic_definition = dict(
        mosaicjson="0.0.1",
        minzoom=minzoom,
        maxzoom=maxzoom,
        bounds=bounds,
        center=[(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom],
        tiles={},
    )

    for tile in tiles:
        z, x, y = list(map(int, tile.split("-")))
        tile = mercantile.Tile(x=x, y=y, z=z)
        quadkey = mercantile.quadkey(*tile)
        geometry = box(*mercantile.bounds(tile))
        intersect_dataset = list(
            filter(lambda x: geometry.intersects(shape(x["geometry"])), dataset)
        )
        if len(intersect_dataset):
            # We limit the item per quadkey to 20
            intersect_dataset = intersect_dataset[0:maximum_items_per_tile]

            mosaic_definition["tiles"][quadkey] = [
                scene["properties"]["landsat:product_id"] for scene in intersect_dataset
            ]

    return mosaic_definition