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

# 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 [3]:

from typing import 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 BaseBackend
from cogeo_mosaic.backends.stac import _fetch, default_stac_accessor
from cogeo_mosaic.mosaic import MosaicJSON


@attr.s
class DynamicStacBackend(BaseBackend):
    """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)

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

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

    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"

    @minzoom.default
    def _minzoom(self):
        return self.tms.minzoom

    @maxzoom.default
    def _maxzoom(self):
        return self.tms.maxzoom

    def __attrs_post_init__(self):
        """Post Init."""
        # 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) -> 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 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 BaseBackend from cogeo_mosaic.backends.stac import _fetch, default_stac_accessor from cogeo_mosaic.mosaic import MosaicJSON @attr.s class DynamicStacBackend(BaseBackend): """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) minzoom: int = attr.ib() maxzoom: int = attr.ib() tms: morecantile.TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS) 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" @minzoom.default def _minzoom(self): return self.tms.minzoom @maxzoom.default def _maxzoom(self): return self.tms.maxzoom def __attrs_post_init__(self): """Post Init.""" # 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) -> 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 [4]:

##  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 [12]:

# 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) center=(0.0, 0.0, 8) minzoom=8 maxzoom=12 name="it's fake but it's ok" quadkeys=[]
5
ImageData(array=masked_array(
  data=[[[283, 350, 339, ..., 536, 439, 386],
         [317, 315, 320, ..., 334, 311, 406],
         [362, 393, 250, ..., 402, 303, 419],
         ...,
         [711, 618, 495, ..., 504, 313, 373],
         [746, 666, 407, ..., 532, 372, 381],
         [273, 569, 506, ..., 546, 606, 418]]],
  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=999999,
  dtype=uint16), assets=['https://earth-search.aws.element84.com/v0/collections/sentinel-s2-l2a-cogs/items/S2A_32UMD_20230419_1_L2A'], bounds=BoundingBox(left=900122.4450862072, bottom=6887893.492833819, right=939258.2035682201, top=6927029.25131583), crs=CRS.from_epsg(3857), metadata={}, band_names=['B01_b1'], dataset_statistics=None, cutline_mask=None)

In [14]:

# 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
S2B_30UXC_20230426_0_L2A [6112]
S2A_30UXC_20230322_0_L2A [442]
S2A_30UXC_20230302_1_L2A [2866]
S2B_30UXC_20230215_0_L2A [409]
S2A_30UXC_20230121_0_L2A [491]

In [ ]: