add ogc maps

pyproject.toml

@@ -32,6 +32,7 @@ dependencies = [
     "cql2>=0.3.6",
     "pydantic>=2.4,<3.0",
     "pydantic-settings~=2.0",
+    "stac-pydantic>=3.0,<4.0",
 ]
 dynamic = ["version"]
 

titiler/pgstac/backend.py

@@ -2,6 +2,9 @@
 
 import json
 import logging
+import math
+import os
+import warnings
 from threading import Lock
 from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Type
 
@@ -15,16 +18,28 @@
 from geojson_pydantic.geometries import Geometry, parse_geometry_obj
 from morecantile import Tile, TileMatrixSet
 from psycopg import errors as pgErrors
+from psycopg.rows import class_row
 from psycopg_pool import ConnectionPool
 from rasterio.crs import CRS
+from rasterio.features import rasterize
 from rasterio.warp import transform, transform_bounds, transform_geom
-from rio_tiler.constants import MAX_THREADS, WEB_MERCATOR_TMS, WGS84_CRS
+from rio_tiler.constants import (
+    MAX_THREADS,
+    WEB_MERCATOR_CRS,
+    WEB_MERCATOR_TMS,
+    WGS84_CRS,
+)
 from rio_tiler.errors import PointOutsideBounds
+from rio_tiler.io import Reader
 from rio_tiler.models import ImageData, PointData
 from rio_tiler.mosaic import mosaic_reader
+
+# _get_width_height, _missing_size were moved in `.utils` in 7.9
+from rio_tiler.reader import _get_width_height, _missing_size
 from rio_tiler.tasks import create_tasks, filter_tasks
 from rio_tiler.types import BBox
 
+from titiler.pgstac.model import Search as TiTilerPGstacSearch
 from titiler.pgstac.reader import SimpleSTACReader
 from titiler.pgstac.settings import CacheSettings, PgstacSettings, RetrySettings
 from titiler.pgstac.utils import retry
@@ -35,6 +50,8 @@
 
 logger = logging.getLogger(__name__)
 
+WORLD_IMG = os.path.join(os.path.dirname(__file__), "data", "world.tif")
+
 
 def multi_points_pgstac(
     asset_list: Sequence[Dict[str, Any]],
@@ -428,3 +445,82 @@ def _reader(item: Dict[str, Any], shape: Dict, **kwargs: Any) -> ImageData:
             **kwargs,
         )
         return img, [x["id"] for x in used_assets]
+
+    def preview(  # noqa: C901
+        self,
+        max_size: Optional[int] = 1024,
+        width: Optional[int] = None,
+        height: Optional[int] = None,
+        dst_crs: Optional[CRS] = None,
+        **kwargs: Any,
+    ) -> Tuple[ImageData, List[str]]:
+        """Create Preview for a Mosaic."""
+        if max_size and (width or height):
+            warnings.warn(
+                "'max_size' will be ignored with with 'height' or 'width' set.",
+                UserWarning,
+                stacklevel=2,
+            )
+            max_size = None
+
+        with self.pool.connection() as conn:
+            with conn.cursor(row_factory=class_row(TiTilerPGstacSearch)) as cursor:
+                cursor.execute(
+                    "SELECT * FROM searches WHERE hash=%s;",
+                    (self.input,),
+                )
+                search_info = cursor.fetchone()
+
+        shapes = [Polygon.from_bounds(-180, -90, 180, 90).model_dump(exclude_none=True)]
+
+        if search_info.metadata.extent and (
+            extent := search_info.metadata.extent.spatial
+        ):
+            shapes = [
+                Polygon.from_bounds(*bbox).model_dump(exclude_none=True)
+                for bbox in extent.bbox
+            ]
+        elif search_info.metadata.bounds:
+            shapes = [
+                Polygon.from_bounds(*search_info.metadata.bounds).model_dump(
+                    exclude_none=True
+                )
+            ]
+
+        with Reader(WORLD_IMG) as src:
+            image = src.read()
+
+        arr = rasterize(
+            shapes,
+            out_shape=(image.height, image.width),
+            transform=image.transform,
+            all_touched=True,
+            default_value=1,
+            fill=0,
+            dtype="uint8",
+        )
+        if not arr.all():
+            image.array[:, arr != 0] = 0
+
+        if dst_crs:
+            src_bounds = list(image.bounds)
+            if dst_crs == WEB_MERCATOR_CRS:
+                src_bounds[1] = max(src_bounds[1], -85.06)
+                src_bounds[3] = min(src_bounds[3], 85.06)
+                image = image.clip(src_bounds)
+            image = image.reproject(dst_crs)
+
+        if max_size:
+            height, width = _get_width_height(max_size, image.height, image.width)
+
+        elif _missing_size(height, width):
+            ratio = image.height / image.width
+            if width:
+                height = math.ceil(width * ratio)
+            else:
+                width = math.ceil(height / ratio)
+
+        if (height and width) and (height != image.height or width != image.width):
+            image = image.resize(height, width, resampling_method="nearest")
+
+        return image, []

titiler/pgstac/data/README.md

@@ -0,0 +1 @@
+The `world.tif` was created using https://www.naturalearthdata.com/downloads/10m-raster-data/10m-natural-earth-2/

titiler/pgstac/dependencies.py

@@ -135,6 +135,7 @@ def get_collection_id(  # noqa: C901
             metadata = model.Metadata(
                 name=f"Mosaic for '{collection_id}' Collection",
                 bounds=search.bbox or collection_bbox[0],
+                extent=collection["extent"],
             )
 
             # item-assets https://github.com/stac-extensions/item-assets

titiler/pgstac/factory.py

@@ -32,7 +32,7 @@
 from psycopg.rows import class_row, dict_row
 from pydantic import Field
 from rio_tiler.constants import MAX_THREADS, WGS84_CRS
-from rio_tiler.utils import CRS_to_urn
+from rio_tiler.utils import CRS_to_uri, CRS_to_urn
 from starlette.datastructures import QueryParams
 from starlette.requests import Request
 from starlette.responses import Response
@@ -45,6 +45,7 @@
     DefaultDependency,
     DstCRSParams,
     HistogramParams,
+    OGCMapsParams,
     PartFeatureParams,
     StatisticsParams,
 )
@@ -113,6 +114,7 @@ class MosaicTilerFactory(BaseFactory):
     add_viewer: bool = False
     add_statistics: bool = False
     add_part: bool = False
+    add_ogc_maps: bool = False
 
     conforms_to: Set[str] = field(
         factory=lambda: {
@@ -142,6 +144,9 @@ def register_routes(self) -> None:
         if self.add_part:
             self.part()
 
+        if self.add_ogc_maps:
+            self.ogc_maps()
+
         if self.add_statistics:
             self.statistics()
 
@@ -694,6 +699,107 @@ def feature_image(
 
             return Response(content, media_type=media_type, headers=headers)
 
+    ############################################################################
+    # OGC Maps (Optional)
+    ############################################################################
+    def ogc_maps(self):  # noqa: C901
+        """Register OGC Maps /map` endpoint."""
+
+        self.conforms_to.update(
+            {
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/core",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/crs",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/scaling",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/scaling/width-definition",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/scaling/height-definition",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/spatial-subsetting",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/spatial-subsetting/bbox-definition",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/spatial-subsetting/bbox-crs",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/spatial-subsetting/crs-curie",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/png",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/jpeg",
+                "https://www.opengis.net/spec/ogcapi-maps-1/1.0/conf/tiff",
+            }
+        )
+
+        # GET endpoints
+        @self.router.get(
+            "/map",
+            operation_id=f"{self.operation_prefix}getMap",
+            **img_endpoint_params,
+        )
+        def get_map(
+            request: Request,
+            search_id=Depends(self.path_dependency),
+            backend_params=Depends(self.backend_dependency),
+            assets_accessor_params=Depends(self.assets_accessor_dependency),
+            reader_params=Depends(self.reader_dependency),
+            ogc_params=Depends(OGCMapsParams),
+            layer_params=Depends(self.layer_dependency),
+            dataset_params=Depends(self.dataset_dependency),
+            pixel_selection=Depends(self.pixel_selection_dependency),
+            post_process=Depends(self.process_dependency),
+            colormap=Depends(self.colormap_dependency),
+            render_params=Depends(self.render_dependency),
+            env=Depends(self.environment_dependency),
+        ):
+            """OGC Maps API."""
+            with rasterio.Env(**env):
+                logger.info(
+                    f"opening data with backend: {self.backend} and reader {self.dataset_reader}"
+                )
+                with self.backend(
+                    search_id,
+                    reader=self.dataset_reader,
+                    reader_options=reader_params.as_dict(),
+                    **backend_params.as_dict(),
+                ) as src_dst:
+                    if ogc_params.bbox is not None:
+                        image, assets = src_dst.part(
+                            ogc_params.bbox,
+                            dst_crs=ogc_params.crs or src_dst.crs,
+                            bounds_crs=ogc_params.bbox_crs or WGS84_CRS,
+                            pixel_selection=pixel_selection,
+                            width=ogc_params.width,
+                            height=ogc_params.height,
+                            max_size=ogc_params.max_size,
+                            **layer_params.as_dict(),
+                            **dataset_params.as_dict(),
+                            **assets_accessor_params.as_dict(),
+                        )
+
+                    else:
+                        image, assets = src_dst.preview(
+                            width=ogc_params.width,
+                            height=ogc_params.height,
+                            max_size=ogc_params.max_size,
+                            dst_crs=ogc_params.crs,
+                            **assets_accessor_params.as_dict(),
+                        )
+
+                    dst_colormap = getattr(src_dst, "colormap", None)
+
+            if post_process:
+                logger.info("post processing image")
+                image = post_process(image)
+
+            content, media_type = self.render_func(
+                image,
+                output_format=ogc_params.format,
+                colormap=colormap or dst_colormap,
+                **render_params.as_dict(),
+            )
+
+            headers: Dict[str, str] = {}
+            if image.bounds is not None:
+                headers["Content-Bbox"] = ",".join(map(str, image.bounds))
+            if uri := CRS_to_uri(image.crs):
+                headers["Content-Crs"] = f"<{uri}>"
+            if OptionalHeader.x_assets in self.optional_headers:
+                headers["X-Assets"] = ",".join(assets)
+
+            return Response(content, media_type=media_type, headers=headers)
+
 
 def add_search_register_route(  # noqa: C901
     app: FastAPI,

titiler/pgstac/main.py

@@ -213,6 +213,7 @@ def pgstac_info(request: Request) -> Dict:
     add_statistics=True,
     add_viewer=True,
     add_part=True,
+    add_ogc_maps=True,
     extensions=[
         searchInfoExtension(),
     ],
@@ -251,6 +252,7 @@ def pgstac_info(request: Request) -> Dict:
     add_statistics=True,
     add_viewer=True,
     add_part=True,
+    add_ogc_maps=True,
     extensions=[
         searchInfoExtension(),
     ],
@@ -268,6 +270,7 @@ def pgstac_info(request: Request) -> Dict:
     path_dependency=ItemIdParams,
     router_prefix="/collections/{collection_id}/items/{item_id}",
     add_viewer=True,
+    add_ogc_maps=True,
     templates=templates,
 )
 app.include_router(

titiler/pgstac/model.py

@@ -13,6 +13,7 @@
 from geojson_pydantic.geometries import Geometry
 from geojson_pydantic.types import BBox
 from pydantic import BaseModel, Field, ValidationInfo, field_validator, model_validator
+from stac_pydantic.collection import Extent
 from typing_extensions import Annotated
 
 from titiler.core.resources.enums import MediaType
@@ -29,6 +30,7 @@ class Metadata(BaseModel):
 
     # WGS84 bounds
     bounds: Optional[BBox] = None
+    extent: Optional[Extent] = None
 
     # Min/Max zoom for WebMercatorQuad TMS
     minzoom: Optional[int] = None