Create a Dynamic STAC backend¶

By default cogeo-mosaic backends were meant to handle writing and reading mosaicjson either from a file or from a database.

While this is fine for most use cases, some users could want something more dynamic. In this Notebook we will show how to create a Dynamic mosaic backend based on STAC api.

Requirements¶

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

• cogeo-mosaic

In [1]:

# Uncomment this line if you need to install the dependencies
# !pip install cogeo-mosaic

# Uncomment this line if you need to install the dependencies # !pip install cogeo-mosaic

In [1]:

from typing import Any, Dict, Tuple, Type, Optional, List

import attr
import morecantile
from rasterio.crs import CRS
from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
from rio_tiler.io import BaseReader
from rio_tiler.io import STACReader

from cogeo_mosaic.backends.base import MosaicJSONBackend
from cogeo_mosaic.backends.stac import _fetch, default_stac_accessor
from cogeo_mosaic.mosaic import MosaicJSON


@attr.s
class DynamicStacBackend(MosaicJSONBackend):
    """Like a STAC backend but dynamic"""

    # input should be the STAC-API url
    input: str = attr.ib()

    # Addition required attribute (STAC Query)
    query: Dict = attr.ib(factory=dict)

    tms: morecantile.TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)

    minzoom: int = attr.ib(default=None)
    maxzoom: int = attr.ib(default=None)

    reader: Type[BaseReader] = attr.ib(default=STACReader)
    reader_options: Dict = attr.ib(factory=dict)

    bounds: Tuple[float, float, float, float] = attr.ib(
        default=(-180, -90, 180, 90)
    )
    crs: CRS = attr.ib(default=WGS84_CRS)

    # STAC API related options
    # max_items |  next_link_key | limit
    stac_api_options: Dict = attr.ib(factory=dict)

    # The reader is read-only, we can't pass mosaic_def to the init method
    mosaic_def: MosaicJSON = attr.ib(init=False)

    _backend_name = "DynamicSTAC"

    def __attrs_post_init__(self):
        """Post Init."""
        self.minzoom = self.minzoom if self.minzoom is not None else self.tms.minzoom
        self.maxzoom = self.maxzoom if self.maxzoom is not None else self.tms.maxzoom
        
        # Construct a FAKE/Empty mosaicJSON
        # mosaic_def has to be defined. As we do for the DynamoDB and SQLite backend
        self.mosaic_def = MosaicJSON(
            mosaicjson="0.0.3",
            name="it's fake but it's ok",
            bounds=self.bounds,
            minzoom=self.minzoom,
            maxzoom=self.maxzoom,
            tiles={}  # we set `tiles` to an empty list.
        )

    def write(self, overwrite: bool = True):
        """This method is not used but is required by the abstract class."""
        raise NotImplementedError

    def update(self):
        """We overwrite the default method."""
        raise NotImplementedError

    def _read(self) -> MosaicJSON:
        """This method is not used but is required by the abstract class."""
        pass

    def assets_for_tile(self, x: int, y: int, z: int) -> List[str]:
        """Retrieve assets for tile."""
        bounds = self.tms.bounds(x, y, z)
        geom = {
            "type": "Polygon",
            "coordinates": [
                [
                    [bounds[0], bounds[3]],
                    [bounds[0], bounds[1]],
                    [bounds[2], bounds[1]],
                    [bounds[2], bounds[3]],
                    [bounds[0], bounds[3]],
                ]
            ],
        }
        return self.get_assets(geom)

    def assets_for_point(self, lng: float, lat: float, **kwargs: Any) -> List[str]:
        """Retrieve assets for point.

        Note: some API only accept Polygon.
        """
        EPSILON = 1e-14
        geom = {
            "type": "Polygon",
            "coordinates": [
                [
                    [lng - EPSILON, lat + EPSILON],
                    [lng - EPSILON, lat - EPSILON],
                    [lng + EPSILON, lat - EPSILON],
                    [lng + EPSILON, lat + EPSILON],
                    [lng - EPSILON, lat + EPSILON],
                ]
            ],
        }
        return self.get_assets(geom)

    def get_assets(self, geom) -> List[str]:
        """Send query to the STAC-API and retrieve assets."""
        query = self.query.copy()
        query["intersects"] = geom

        features = _fetch(
            self.input,
            query,
            **self.stac_api_options,
        )
        return [default_stac_accessor(f) for f in features]

    @property
    def _quadkeys(self) -> List[str]:
        return []

from typing import Any, Dict, Tuple, Type, Optional, List import attr import morecantile from rasterio.crs import CRS from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS from rio_tiler.io import BaseReader from rio_tiler.io import STACReader from cogeo_mosaic.backends.base import MosaicJSONBackend from cogeo_mosaic.backends.stac import _fetch, default_stac_accessor from cogeo_mosaic.mosaic import MosaicJSON @attr.s class DynamicStacBackend(MosaicJSONBackend): """Like a STAC backend but dynamic""" # input should be the STAC-API url input: str = attr.ib() # Addition required attribute (STAC Query) query: Dict = attr.ib(factory=dict) tms: morecantile.TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS) minzoom: int = attr.ib(default=None) maxzoom: int = attr.ib(default=None) reader: Type[BaseReader] = attr.ib(default=STACReader) reader_options: Dict = attr.ib(factory=dict) bounds: Tuple[float, float, float, float] = attr.ib( default=(-180, -90, 180, 90) ) crs: CRS = attr.ib(default=WGS84_CRS) # STAC API related options # max_items | next_link_key | limit stac_api_options: Dict = attr.ib(factory=dict) # The reader is read-only, we can't pass mosaic_def to the init method mosaic_def: MosaicJSON = attr.ib(init=False) _backend_name = "DynamicSTAC" def __attrs_post_init__(self): """Post Init.""" self.minzoom = self.minzoom if self.minzoom is not None else self.tms.minzoom self.maxzoom = self.maxzoom if self.maxzoom is not None else self.tms.maxzoom # Construct a FAKE/Empty mosaicJSON # mosaic_def has to be defined. As we do for the DynamoDB and SQLite backend self.mosaic_def = MosaicJSON( mosaicjson="0.0.3", name="it's fake but it's ok", bounds=self.bounds, minzoom=self.minzoom, maxzoom=self.maxzoom, tiles={} # we set `tiles` to an empty list. ) def write(self, overwrite: bool = True): """This method is not used but is required by the abstract class.""" raise NotImplementedError def update(self): """We overwrite the default method.""" raise NotImplementedError def _read(self) -> MosaicJSON: """This method is not used but is required by the abstract class.""" pass def assets_for_tile(self, x: int, y: int, z: int) -> List[str]: """Retrieve assets for tile.""" bounds = self.tms.bounds(x, y, z) geom = { "type": "Polygon", "coordinates": [ [ [bounds[0], bounds[3]], [bounds[0], bounds[1]], [bounds[2], bounds[1]], [bounds[2], bounds[3]], [bounds[0], bounds[3]], ] ], } return self.get_assets(geom) def assets_for_point(self, lng: float, lat: float, **kwargs: Any) -> List[str]: """Retrieve assets for point. Note: some API only accept Polygon. """ EPSILON = 1e-14 geom = { "type": "Polygon", "coordinates": [ [ [lng - EPSILON, lat + EPSILON], [lng - EPSILON, lat - EPSILON], [lng + EPSILON, lat - EPSILON], [lng + EPSILON, lat + EPSILON], [lng - EPSILON, lat + EPSILON], ] ], } return self.get_assets(geom) def get_assets(self, geom) -> List[str]: """Send query to the STAC-API and retrieve assets.""" query = self.query.copy() query["intersects"] = geom features = _fetch( self.input, query, **self.stac_api_options, ) return [default_stac_accessor(f) for f in features] @property def _quadkeys(self) -> List[str]: return []

In [2]:

##  Base Query for sat-api
# - limit of 5 items per page (we will stop at page 1)
# - less than 25% of cloud
# - more than 75% of data coverage
# - `sentinel-s2-l2a-cogs` collection
query = {
    "collections": ["sentinel-s2-l2a-cogs"],
    "limit": 5,
    "query": {
        "eo:cloud_cover": {
            "lt": 25
        },
        "sentinel:data_coverage": {
            "gt": 75
        }
    },
    "fields": {
      'include': ['id'],
      'exclude': ['assets', 'geometry']
    }
}

## Base Query for sat-api # - limit of 5 items per page (we will stop at page 1) # - less than 25% of cloud # - more than 75% of data coverage # - `sentinel-s2-l2a-cogs` collection query = { "collections": ["sentinel-s2-l2a-cogs"], "limit": 5, "query": { "eo:cloud_cover": { "lt": 25 }, "sentinel:data_coverage": { "gt": 75 } }, "fields": { 'include': ['id'], 'exclude': ['assets', 'geometry'] } }

In [3]:

# Read Tile
with DynamicStacBackend(
    "https://earth-search.aws.element84.com/v0/search",
    query=query,
    stac_api_options={"max_items": 5},
    minzoom=8,  # we know this by analysing the sentinel 2 data
    maxzoom=12,  # we know this by analysing the sentinel 2 data
) as mosaic:
    print(mosaic.info())
    print(len(mosaic.assets_for_tile(535, 335, 10)))
    img, _ = mosaic.tile(535, 335, 10, assets="B01", tilesize=128, threads=0)

print(img)

# Read Tile with DynamicStacBackend( "https://earth-search.aws.element84.com/v0/search", query=query, stac_api_options={"max_items": 5}, minzoom=8, # we know this by analysing the sentinel 2 data maxzoom=12, # we know this by analysing the sentinel 2 data ) as mosaic: print(mosaic.info()) print(len(mosaic.assets_for_tile(535, 335, 10))) img, _ = mosaic.tile(535, 335, 10, assets="B01", tilesize=128, threads=0) print(img)

bounds=(-180.0, -90.0, 180.0, 90.0) crs='http://www.opengis.net/def/crs/EPSG/0/4326' name="it's fake but it's ok" quadkeys=[] mosaic_tilematrixset="<TileMatrixSet title='Google Maps Compatible for the World' id='WebMercatorQuad' crs='http://www.opengis.net/def/crs/EPSG/0/3857>" mosaic_minzoom=8 mosaic_maxzoom=12
5
ImageData(array=masked_array(
  data=[[[3964, 3867, 4477, ..., 7267, 6752, 6690],
         [4358, 4450, 4357, ..., 5665, 5759, 5270],
         [4624, 5443, 5051, ..., 7196, 7560, 6267],
         ...,
         [6848, 5832, 5926, ..., 6824, 7578, 5917],
         [5687, 5581, 5709, ..., 6956, 6578, 6183],
         [5896, 5798, 5450, ..., 6242, 6246, 6814]]],
  mask=[[[False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         ...,
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False],
         [False, False, False, ..., False, False, False]]],
  fill_value=np.uint64(999999),
  dtype=uint16), assets=['https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/S2C_32UMD_20251107_0_L2A'], bounds=BoundingBox(left=900122.445086237, bottom=6887893.4928338025, right=939258.2035682462, top=6927029.2513158135), crs=CRS.from_wkt('PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]'), metadata={'mosaic_method': 'FirstMethod', 'mosaic_assets_count': 5, 'mosaic_assets_used': 1, 'timings': [('search', 0.51), ('mosaicking', 3405.62)]}, nodata=None, scales=[1.0], offsets=[0.0], band_names=['B01_b1'], band_descriptions=['B01_'], dataset_statistics=None, cutline_mask=None, alpha_mask=None)

In [23]:

# Read Point values
with DynamicStacBackend(
    "https://earth-search.aws.element84.com/v0/search",
    query=query,
    stac_api_options={"max_items": 5},
    minzoom=8,  # we know this by analysing the sentinel 2 data
    maxzoom=12,  # we know this by analysing the sentinel 2 data
) as mosaic:
    print(len(mosaic.assets_for_point(-1.0546875,  51.99)))
    values = mosaic.point(-1.0546875,  51.99, assets="B01")

    for (f, v) in values:
        print(f.split("/")[-1], v.data)

# Read Point values with DynamicStacBackend( "https://earth-search.aws.element84.com/v0/search", query=query, stac_api_options={"max_items": 5}, minzoom=8, # we know this by analysing the sentinel 2 data maxzoom=12, # we know this by analysing the sentinel 2 data ) as mosaic: print(len(mosaic.assets_for_point(-1.0546875, 51.99))) values = mosaic.point(-1.0546875, 51.99, assets="B01") for (f, v) in values: print(f.split("/")[-1], v.data)

5
S2A_30UXC_20251118_0_L2A [241]
S2B_30UXC_20251101_0_L2A [5133]
S2B_30UXC_20250922_1_L2A [318]
S2A_30UXC_20250711_1_L2A [525]
S2C_30UXC_20250619_0_L2A [362]