add support for the STAC projection extension

CHANGES.md

@@ -4,6 +4,10 @@
 * Enable dynamic definition of Asset **reader** in `MultiBaseReader` (https://github.com/cogeotiff/rio-tiler/pull/711/, https://github.com/cogeotiff/rio-tiler/pull/728)
 * Adding `default_assets` for MultiBaseReader and STACReader (author @mccarthyryanc, https://github.com/cogeotiff/rio-tiler/pull/722)
 * Adding `default_bands` for MultiBandReader (https://github.com/cogeotiff/rio-tiler/pull/722)
+* Adding support for the STAC `Projection` extension to derive the `bounds`, `crs`, `minzoom` and `maxzoom` properties  **breaking change**
+* Refactor internal function and base classes for the `minzoom/maxzoom` calculation **breaking change**
+* Adding `transform`, `height` and `width` attributes (outside init) for `SpatialMixin` class
+* Moved `_dst_geom_in_tms_crs` from Reader to `SpatialMixin` class **breaking change**
 
 # 6.7.0 (2024-09-05)
 

rio_tiler/io/base.py

@@ -9,9 +9,11 @@
 
 import attr
 import numpy
+from affine import Affine
 from morecantile import Tile, TileMatrixSet
 from rasterio.crs import CRS
-from rasterio.warp import transform_bounds
+from rasterio.rio.overview import get_maximum_overview_level
+from rasterio.warp import calculate_default_transform, transform_bounds
 
 from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
 from rio_tiler.errors import (
@@ -42,6 +44,10 @@ class SpatialMixin:
     bounds: BBox = attr.ib(init=False)
     crs: CRS = attr.ib(init=False)
 
+    transform: Optional[Affine] = attr.ib(default=None, init=False)
+    height: Optional[int] = attr.ib(default=None, init=False)
+    width: Optional[int] = attr.ib(default=None, init=False)
+
     geographic_crs: CRS = attr.ib(init=False, default=WGS84_CRS)
 
     @cached_property
@@ -85,6 +91,77 @@ def geographic_bounds(self) -> BBox:
 
         return bounds
 
+    @cached_property
+    def _dst_geom_in_tms_crs(self):
+        """Return dataset geom info in TMS projection."""
+        tms_crs = self.tms.rasterio_crs
+        if self.crs != tms_crs:
+            dst_affine, w, h = calculate_default_transform(
+                self.crs,
+                tms_crs,
+                self.width,
+                self.height,
+                *self.bounds,
+            )
+        else:
+            dst_affine = list(self.transform)
+            w = self.width
+            h = self.height
+
+        return dst_affine, w, h
+
+    @cached_property
+    def _minzoom(self) -> int:
+        """Calculate dataset minimum zoom level."""
+        # We assume the TMS tilesize to be constant over all matrices
+        # ref: https://github.com/OSGeo/gdal/blob/dc38aa64d779ecc45e3cd15b1817b83216cf96b8/gdal/frmts/gtiff/cogdriver.cpp#L274
+        tilesize = self.tms.tileMatrices[0].tileWidth
+
+        if all([self.transform, self.height, self.width]):
+            try:
+                dst_affine, w, h = self._dst_geom_in_tms_crs
+
+                # The minzoom is defined by the resolution of the maximum theoretical overview level
+                # We assume `tilesize`` is the smallest overview size
+                overview_level = get_maximum_overview_level(w, h, minsize=tilesize)
+
+                # Get the resolution of the overview
+                resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
+                ovr_resolution = resolution * (2**overview_level)
+
+                # Find what TMS matrix match the overview resolution
+                return self.tms.zoom_for_res(ovr_resolution)
+
+            except:  # noqa
+                # if we can't get max zoom from the dataset we default to TMS maxzoom
+                warnings.warn(
+                    "Cannot determine minzoom based on dataset information, will default to TMS minzoom.",
+                    UserWarning,
+                )
+
+        return self.tms.minzoom
+
+    @cached_property
+    def _maxzoom(self) -> int:
+        """Calculate dataset maximum zoom level."""
+        if all([self.transform, self.height, self.width]):
+            try:
+                dst_affine, _, _ = self._dst_geom_in_tms_crs
+
+                # The maxzoom is defined by finding the minimum difference between
+                # the raster resolution and the zoom level resolution
+                resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
+                return self.tms.zoom_for_res(resolution)
+
+            except:  # noqa
+                # if we can't get min/max zoom from the dataset we default to TMS maxzoom
+                warnings.warn(
+                    "Cannot determine maxzoom based on dataset information, will default to TMS maxzoom.",
+                    UserWarning,
+                )
+
+        return self.tms.maxzoom
+
     def tile_exists(self, tile_x: int, tile_y: int, tile_z: int) -> bool:
         """Check if a tile intersects the dataset bounds.
 
@@ -267,7 +344,6 @@ class MultiBaseReader(SpatialMixin, metaclass=abc.ABCMeta):
     reader_options: Dict = attr.ib(factory=dict)
 
     assets: Sequence[str] = attr.ib(init=False)
-
     default_assets: Optional[Sequence[str]] = attr.ib(init=False, default=None)
 
     ctx: Any = attr.ib(init=False, default=contextlib.nullcontext)
@@ -945,7 +1021,6 @@ class MultiBandReader(SpatialMixin, metaclass=abc.ABCMeta):
     reader_options: Dict = attr.ib(factory=dict)
 
     bands: Sequence[str] = attr.ib(init=False)
-
     default_bands: Optional[Sequence[str]] = attr.ib(init=False, default=None)
 
     def __enter__(self):

rio_tiler/io/rasterio.py

@@ -19,7 +19,7 @@
 from rasterio.rio.overview import get_maximum_overview_level
 from rasterio.transform import from_bounds as transform_from_bounds
 from rasterio.vrt import WarpedVRT
-from rasterio.warp import calculate_default_transform, transform_geom
+from rasterio.warp import transform_geom
 from rasterio.windows import Window
 from rasterio.windows import from_bounds as window_from_bounds
 
@@ -83,8 +83,6 @@ class Reader(BaseReader):
 
     # Context Manager to handle rasterio open/close
     _ctx_stack = attr.ib(init=False, factory=contextlib.ExitStack)
-    _minzoom: int = attr.ib(init=False, default=None)
-    _maxzoom: int = attr.ib(init=False, default=None)
 
     @options.default
     def _options_default(self):
@@ -121,6 +119,10 @@ def __attrs_post_init__(self):
         self.bounds = tuple(self.dataset.bounds)
         self.crs = self.dataset.crs
 
+        self.transform = self.dataset.transform
+        self.height = self.dataset.height
+        self.width = self.dataset.width
+
         if self.colormap is None:
             self._get_colormap()
 
@@ -140,85 +142,15 @@ def __exit__(self, exc_type, exc_value, traceback):
         """Support using with Context Managers."""
         self.close()
 
-    def _dst_geom_in_tms_crs(self):
-        """Return dataset info in TMS projection."""
-        tms_crs = self.tms.rasterio_crs
-        if self.crs != tms_crs:
-            dst_affine, w, h = calculate_default_transform(
-                self.crs,
-                tms_crs,
-                self.dataset.width,
-                self.dataset.height,
-                *self.dataset.bounds,
-            )
-        else:
-            dst_affine = list(self.dataset.transform)
-            w = self.dataset.width
-            h = self.dataset.height
-
-        return dst_affine, w, h
-
-    def get_minzoom(self) -> int:
-        """Define dataset minimum zoom level."""
-        if self._minzoom is None:
-            # We assume the TMS tilesize to be constant over all matrices
-            # ref: https://github.com/OSGeo/gdal/blob/dc38aa64d779ecc45e3cd15b1817b83216cf96b8/gdal/frmts/gtiff/cogdriver.cpp#L274
-            tilesize = self.tms.tileMatrices[0].tileWidth
-
-            try:
-                dst_affine, w, h = self._dst_geom_in_tms_crs()
-
-                # The minzoom is defined by the resolution of the maximum theoretical overview level
-                # We assume `tilesize`` is the smallest overview size
-                overview_level = get_maximum_overview_level(w, h, minsize=tilesize)
-
-                # Get the resolution of the overview
-                resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
-                ovr_resolution = resolution * (2**overview_level)
-
-                # Find what TMS matrix match the overview resolution
-                self._minzoom = self.tms.zoom_for_res(ovr_resolution)
-
-            except:  # noqa
-                # if we can't get max zoom from the dataset we default to TMS maxzoom
-                warnings.warn(
-                    "Cannot determine minzoom based on dataset information, will default to TMS minzoom.",
-                    UserWarning,
-                )
-                self._minzoom = self.tms.minzoom
-
-        return self._minzoom
-
-    def get_maxzoom(self) -> int:
-        """Define dataset maximum zoom level."""
-        if self._maxzoom is None:
-            try:
-                dst_affine, _, _ = self._dst_geom_in_tms_crs()
-
-                # The maxzoom is defined by finding the minimum difference between
-                # the raster resolution and the zoom level resolution
-                resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
-                self._maxzoom = self.tms.zoom_for_res(resolution)
-
-            except:  # noqa
-                # if we can't get min/max zoom from the dataset we default to TMS maxzoom
-                warnings.warn(
-                    "Cannot determine maxzoom based on dataset information, will default to TMS maxzoom.",
-                    UserWarning,
-                )
-                self._maxzoom = self.tms.maxzoom
-
-        return self._maxzoom
-
     @property
     def minzoom(self):
         """Return dataset minzoom."""
-        return self.get_minzoom()
+        return self._minzoom
 
     @property
     def maxzoom(self):
         """Return dataset maxzoom."""
-        return self.get_maxzoom()
+        return self._maxzoom
 
     def _get_colormap(self):
         """Retrieve the internal colormap."""
@@ -676,23 +608,30 @@ def __attrs_post_init__(self):
 
         height, width = self.dataset.height, self.dataset.width
         self.bounds = (0, height, width, 0)
-        self.transform = transform_from_bounds(*self.bounds, width=width, height=height)
-
         self.tms = LocalTileMatrixSet(width=width, height=height)
-        self._minzoom = self.tms.minzoom
-        self._maxzoom = self.tms.maxzoom
+        self.transform = transform_from_bounds(*self.bounds, width=width, height=height)
+        self.height = height
+        self.width = width
 
         if self.colormap is None:
             self._get_colormap()
 
-        if min(
-            self.dataset.width, self.dataset.height
-        ) > 512 and not self.dataset.overviews(1):
+        if min(width, height) > 512 and not self.dataset.overviews(1):
             warnings.warn(
                 "The dataset has no Overviews. rio-tiler performances might be impacted.",
                 NoOverviewWarning,
             )
 
+    @property
+    def minzoom(self):
+        """Return dataset minzoom."""
+        return self.tms.minzoom
+
+    @property
+    def maxzoom(self):
+        """Return dataset maxzoom."""
+        return self.tms.maxzoom
+
     def tile(  # type: ignore
         self,
         tile_x: int,

rio_tiler/io/stac.py

@@ -14,6 +14,7 @@
 from cachetools.keys import hashkey
 from morecantile import TileMatrixSet
 from rasterio.crs import CRS
+from rasterio.transform import array_bounds
 
 from rio_tiler.constants import WEB_MERCATOR_TMS, WGS84_CRS
 from rio_tiler.errors import InvalidAssetName, MissingAssets
@@ -228,8 +229,8 @@ class STACReader(MultiBaseReader):
     item: pystac.Item = attr.ib(default=None, converter=_to_pystac_item)
 
     tms: TileMatrixSet = attr.ib(default=WEB_MERCATOR_TMS)
-    minzoom: int = attr.ib()
-    maxzoom: int = attr.ib()
+    minzoom: int = attr.ib(default=None)
+    maxzoom: int = attr.ib(default=None)
 
     geographic_crs: CRS = attr.ib(default=WGS84_CRS)
 
@@ -239,6 +240,7 @@ class STACReader(MultiBaseReader):
     include_asset_types: Set[str] = attr.ib(default=DEFAULT_VALID_TYPE)
     exclude_asset_types: Optional[Set[str]] = attr.ib(default=None)
 
+    assets: Sequence[str] = attr.ib(init=False)
     default_assets: Optional[Sequence[str]] = attr.ib(default=None)
 
     reader: Type[BaseReader] = attr.ib(default=Reader)
@@ -254,10 +256,25 @@ def __attrs_post_init__(self):
             fetch(self.input, **self.fetch_options), self.input
         )
 
-        # TODO: get bounds/crs using PROJ extension if available
         self.bounds = self.item.bbox
         self.crs = WGS84_CRS
 
+        if self.item.ext.has("proj"):
+            if all(
+                [
+                    self.item.ext.proj.transform,
+                    self.item.ext.proj.shape,
+                    self.item.ext.proj.crs_string,
+                ]
+            ):
+                self.height, self.width = self.item.ext.proj.shape
+                self.transform = self.item.ext.proj.transform
+                self.bounds = array_bounds(self.height, self.width, self.transform)
+                self.crs = rasterio.crs.CRS.from_string(self.item.ext.proj.crs_string)
+
+        self.minzoom = self.minzoom if self.minzoom is not None else self._minzoom
+        self.maxzoom = self.maxzoom if self.maxzoom is not None else self._maxzoom
+
         self.assets = list(
             _get_assets(
                 self.item,
@@ -270,14 +287,6 @@ def __attrs_post_init__(self):
         if not self.assets:
             raise MissingAssets("No valid asset found. Asset's media types not supported")
 
-    @minzoom.default
-    def _minzoom(self):
-        return self.tms.minzoom
-
-    @maxzoom.default
-    def _maxzoom(self):
-        return self.tms.maxzoom
-
     def _get_reader(self, asset_info: AssetInfo) -> Tuple[Type[BaseReader], Dict]:
         """Get Asset Reader."""
         return self.reader, {}