add coord-crs options (#632)

* add coord-crs options

* fix test

* avoid test to fail for encoding reason

Author: vincentsarago

CHANGES.md

@@ -1,8 +1,9 @@
 # Release Notes
 
-## Next
+## 0.11.6 (2023-04-14)
 
 * Allow a default `rescale` parameter to be set via a dependency (author @samn, https://github.com/developmentseed/titiler/pull/619)
+* add `coord-crs` parameter for `/point`, `/part` and `/feature` endpoints
 
 ## 0.11.5 (2023-03-22)
 

docs/src/endpoints/cog.md

@@ -113,6 +113,7 @@ Example:
     - **url** (str): Cloud Optimized GeoTIFF URL. **Required**
     - **bidx** (array[int]): Dataset band indexes (e.g `bidx=1`, `bidx=1&bidx=2&bidx=3`).
     - **expression** (str): rio-tiler's band math expression (e.g B1/B2).
+    - **coord-crs** (str): Coordinate Reference System of the input coordinates. Default to `epsg:4326`.
     - **max_size** (int): Max image size, default is 1024.
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
@@ -148,6 +149,7 @@ Example:
     - **url** (str): Cloud Optimized GeoTIFF URL. **Required**
     - **bidx** (array[int]): Dataset band indexes (e.g `bidx=1`, `bidx=1&bidx=2&bidx=3`).
     - **expression** (str): rio-tiler's band math expression (e.g B1/B2).
+    - **coord-crs** (str): Coordinate Reference System of the input geometry coordinates. Default to `epsg:4326`.
     - **max_size** (int): Max image size, default is 1024.
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
@@ -182,6 +184,7 @@ Note: if `height` and `width` are provided `max_size` will be ignored.
     - **url** (str): Cloud Optimized GeoTIFF URL. **Required**
     - **bidx** (array[int]): Dataset band indexes (e.g `bidx=1`, `bidx=1&bidx=2&bidx=3`).
     - **expression** (str): rio-tiler's band math expression (e.g B1/B2).
+    - **coord-crs** (str): Coordinate Reference System of the input coordinates. Default to `epsg:4326`.
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
     - **resampling** (str): rasterio resampling method. Default is `nearest`.
@@ -320,6 +323,7 @@ Example:
     - **url** (str): Cloud Optimized GeoTIFF URL. **Required**
     - **bidx** (array[int]): Dataset band indexes (e.g `bidx=1`, `bidx=1&bidx=2&bidx=3`).
     - **expression** (str): rio-tiler's band math expression (e.g B1/B2).
+    - **coord-crs** (str): Coordinate Reference System of the input geometry coordinates. Default to `epsg:4326`.
     - **max_size** (int): Max image size from which to calculate statistics, default is 1024.
     - **height** (int): Force image height from which to calculate statistics.
     - **width** (int): Force image width from which to calculate statistics.

docs/src/endpoints/stac.md

@@ -124,6 +124,7 @@ Example:
     - **expression** (str): rio-tiler's math expression with asset names (e.g `Asset1_b1/Asset2_b1`).
     - **asset_as_band** (bool): tell rio-tiler that each asset is a 1 band dataset, so expression `Asset1/Asset2` can be passed.
     - **asset_bidx** (array[str]): Per asset band math expression (e.g `Asset1|1,2,3`).
+    - **coord-crs** (str): Coordinate Reference System of the input coordinates. Default to `epsg:4326`.
     - **max_size** (int): Max image size, default is 1024.
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
@@ -162,6 +163,7 @@ Example:
     - **expression** (str): rio-tiler's math expression with asset names (e.g `Asset1_b1/Asset2_b1`).
     - **asset_as_band** (bool): tell rio-tiler that each asset is a 1 band dataset, so expression `Asset1/Asset2` can be passed.
     - **asset_bidx** (array[str]): Per asset band math expression (e.g `Asset1|1,2,3`).
+    - **coord-crs** (str): Coordinate Reference System of the input geometry coordinates. Default to `epsg:4326`.
     - **max_size** (int): Max image size, default is 1024.
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
@@ -198,6 +200,7 @@ Example:
     - **expression** (str): rio-tiler's math expression with asset names (e.g `Asset1_b1/Asset2_b1`).
     - **asset_as_band** (bool): tell rio-tiler that each asset is a 1 band dataset, so expression `Asset1/Asset2` can be passed.
     - **asset_bidx** (array[str]): Per asset band math expression (e.g `Asset1|1,2,3`).
+    - **coord-crs** (str): Coordinate Reference System of the input coordinates. Default to `epsg:4326`.
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
     - **resampling** (str): rasterio resampling method. Default is `nearest`.
@@ -389,6 +392,7 @@ Example:
     - **expression** (str): rio-tiler's math expression with asset names (e.g `Asset1_b1/Asset2_b1`).
     - **asset_as_band** (bool): tell rio-tiler that each asset is a 1 band dataset, so expression `Asset1/Asset2` can be passed.
     - **asset_bidx** (array[str]): Per asset band math expression (e.g `Asset1|1,2,3`).
+    - **coord-crs** (str): Coordinate Reference System of the input geometry coordinates. Default to `epsg:4326`.
     - **max_size** (int): Max image size from which to calculate statistics, default is 1024.
     - **height** (int): Force image height from which to calculate statistics.
     - **width** (int): Force image width from which to calculate statistics.

src/titiler/application/tests/routes/test_cog.py

@@ -59,7 +59,7 @@ def test_wmts(rio, app):
     assert response.headers["content-type"] == "application/xml"
     assert response.headers["Cache-Control"] == "private, max-age=3600"
     assert (
-        "http://testserver/cog/WMTSCapabilities.xml?url=https://myurl.com/cog.tif"
+        "http://testserver/cog/WMTSCapabilities.xml?url=https"
         in response.content.decode()
     )
     assert "<ows:Identifier>cogeo</ows:Identifier>" in response.content.decode()

src/titiler/core/tests/test_dependencies.py

@@ -250,10 +250,10 @@ def _assets_bidx(params=Depends(dependencies.AssetsBidxParams)):
     assert response.json()["asset_indexes"] == {"data": [1, 2, 3], "image": [1]}
 
     response = client.get(
-        "/third?assets=data&assets=image&asset_expression=data|b1\b2&asset_expression=image|b1*b2"
+        "/third?assets=data&assets=image&asset_expression=data|b1/b2&asset_expression=image|b1*b2"
     )
     assert response.json()["assets"] == ["data", "image"]
-    assert response.json()["asset_expression"] == {"data": "b1\b2", "image": "b1*b2"}
+    assert response.json()["asset_expression"] == {"data": "b1/b2", "image": "b1*b2"}
 
 
 def test_bands():

src/titiler/core/tests/test_factories.py

@@ -189,6 +189,14 @@ def test_TilerFactory():
     assert len(response.json()["values"]) == 1
     assert response.json()["band_names"] == ["b1"]
 
+    response = client.get(
+        f"/point/-6259272.328324187,12015838.020930404?url={DATA_DIR}/cog.tif&coord-crs=EPSG:3857"
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/json"
+    assert len(response.json()["values"]) == 1
+    assert response.json()["band_names"] == ["b1"]
+
     response = client.get(f"/point/-56.228,72.715?url={DATA_DIR}/cog.tif&bidx=1&bidx=1")
     assert len(response.json()["values"]) == 2
     assert response.json()["band_names"] == ["b1", "b1"]

src/titiler/core/titiler/core/dependencies.py

@@ -7,6 +7,7 @@
 
 import numpy
 from fastapi import HTTPException, Query
+from rasterio.crs import CRS
 from rasterio.enums import Resampling
 from rio_tiler.colormap import cmap, parse_color
 from rio_tiler.errors import MissingAssets, MissingBands
@@ -448,3 +449,17 @@ def __post_init__(self):
 
         if self.range:
             self.range = list(map(float, self.range.split(",")))  # type: ignore
+
+
+def CRSParams(
+    crs: str = Query(
+        None,
+        alias="coord-crs",
+        description="Coordinate Reference System of the input coords. Default to `epsg:4326`.",
+    )
+) -> Optional[CRS]:
+    """Coordinate Reference System Coordinates Param."""
+    if crs:
+        return CRS.from_user_input(crs)
+
+    return None

src/titiler/core/titiler/core/factory.py

@@ -14,6 +14,8 @@
 from morecantile import TileMatrixSet
 from morecantile import tms as morecantile_tms
 from morecantile.defaults import TileMatrixSets
+from rasterio.crs import CRS
+from rio_tiler.constants import WGS84_CRS
 from rio_tiler.io import BaseReader, MultiBandReader, MultiBaseReader, Reader
 from rio_tiler.models import BandStatistics, Bounds, Info
 from rio_tiler.types import ColorMapType
@@ -35,6 +37,7 @@
     BandsParams,
     BidxExprParams,
     ColorMapParams,
+    CRSParams,
     DatasetParams,
     DatasetPathParams,
     DefaultDependency,
@@ -367,6 +370,7 @@ def info_geojson(
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
                     return Feature(
+                        type="Feature",
                         geometry=Polygon.from_bounds(*src_dst.geographic_bounds),
                         properties=src_dst.info(),
                     )
@@ -428,6 +432,7 @@ def geojson_statistics(
                 ..., description="GeoJSON Feature or FeatureCollection."
             ),
             src_path=Depends(self.path_dependency),
+            coord_crs: Optional[CRS] = Depends(CRSParams),
             layer_params=Depends(self.layer_dependency),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_dependency),
@@ -439,13 +444,14 @@ def geojson_statistics(
             """Get Statistics from a geojson feature or featureCollection."""
             fc = geojson
             if isinstance(fc, Feature):
-                fc = FeatureCollection(features=[geojson])
+                fc = FeatureCollection(type="FeatureCollection", features=[geojson])
 
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
                     for feature in fc:
                         data = src_dst.feature(
                             feature.dict(exclude_none=True),
+                            shape_crs=coord_crs or WGS84_CRS,
                             **layer_params,
                             **image_params,
                             **dataset_params,
@@ -845,17 +851,20 @@ def point(
             lon: float = Path(..., description="Longitude"),
             lat: float = Path(..., description="Latitude"),
             src_path=Depends(self.path_dependency),
+            coord_crs: Optional[CRS] = Depends(CRSParams),
             layer_params=Depends(self.layer_dependency),
             dataset_params=Depends(self.dataset_dependency),
             reader_params=Depends(self.reader_dependency),
             env=Depends(self.environment_dependency),
         ):
             """Get Point value for a dataset."""
+
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
                     pts = src_dst.point(
                         lon,
                         lat,
+                        coord_crs=coord_crs or WGS84_CRS,
                         **layer_params,
                         **dataset_params,
                     )
@@ -949,6 +958,7 @@ def part(
             maxy: float = Path(..., description="Bounding box max Y"),
             format: ImageType = Query(..., description="Output image type."),
             src_path=Depends(self.path_dependency),
+            coord_crs: Optional[CRS] = Depends(CRSParams),
             layer_params=Depends(self.layer_dependency),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_dependency),
@@ -969,6 +979,7 @@ def part(
                 with self.reader(src_path, **reader_params) as src_dst:
                     image = src_dst.part(
                         [minx, miny, maxx, maxy],
+                        bounds_crs=coord_crs or WGS84_CRS,
                         **layer_params,
                         **image_params,
                         **dataset_params,
@@ -1014,6 +1025,7 @@ def geojson_crop(
                 None, description="Output image type. Default is auto."
             ),
             src_path=Depends(self.path_dependency),
+            coord_crs: Optional[CRS] = Depends(CRSParams),
             layer_params=Depends(self.layer_dependency),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_dependency),
@@ -1034,6 +1046,7 @@ def geojson_crop(
                 with self.reader(src_path, **reader_params) as src_dst:
                     image = src_dst.feature(
                         geojson.dict(exclude_none=True),
+                        shape_crs=coord_crs or WGS84_CRS,
                         **layer_params,
                         **image_params,
                         **dataset_params,
@@ -1135,6 +1148,7 @@ def info_geojson(
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
                     return Feature(
+                        type="Feature",
                         geometry=Polygon.from_bounds(*src_dst.geographic_bounds),
                         properties={
                             asset: asset_info
@@ -1254,6 +1268,7 @@ def geojson_statistics(
                 ..., description="GeoJSON Feature or FeatureCollection."
             ),
             src_path=Depends(self.path_dependency),
+            coord_crs: Optional[CRS] = Depends(CRSParams),
             layer_params=Depends(AssetsBidxExprParamsOptional),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_dependency),
@@ -1265,7 +1280,7 @@ def geojson_statistics(
             """Get Statistics from a geojson feature or featureCollection."""
             fc = geojson
             if isinstance(fc, Feature):
-                fc = FeatureCollection(features=[geojson])
+                fc = FeatureCollection(type="FeatureCollection", features=[geojson])
 
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
@@ -1276,6 +1291,7 @@ def geojson_statistics(
                     for feature in fc:
                         data = src_dst.feature(
                             feature.dict(exclude_none=True),
+                            shape_crs=coord_crs or WGS84_CRS,
                             **layer_params,
                             **image_params,
                             **dataset_params,
@@ -1370,6 +1386,7 @@ def info_geojson(
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
                     return Feature(
+                        type="Feature",
                         geometry=Polygon.from_bounds(*src_dst.geographic_bounds),
                         properties=src_dst.info(**bands_params),
                     )
@@ -1444,6 +1461,7 @@ def geojson_statistics(
                 ..., description="GeoJSON Feature or FeatureCollection."
             ),
             src_path=Depends(self.path_dependency),
+            coord_crs: Optional[CRS] = Depends(CRSParams),
             bands_params=Depends(BandsExprParamsOptional),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_dependency),
@@ -1455,7 +1473,7 @@ def geojson_statistics(
             """Get Statistics from a geojson feature or featureCollection."""
             fc = geojson
             if isinstance(fc, Feature):
-                fc = FeatureCollection(features=[geojson])
+                fc = FeatureCollection(type="FeatureCollection", features=[geojson])
 
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
@@ -1466,6 +1484,7 @@ def geojson_statistics(
                     for feature in fc:
                         data = src_dst.feature(
                             feature.dict(exclude_none=True),
+                            shape_crs=coord_crs or WGS84_CRS,
                             **bands_params,
                             **image_params,
                             **dataset_params,

src/titiler/mosaic/titiler/mosaic/factory.py

@@ -205,7 +205,9 @@ def info_geojson(
                 ) as src_dst:
                     info = src_dst.info()
                     return Feature(
-                        geometry=Polygon.from_bounds(*info.bounds), properties=info
+                        type="Feature",
+                        geometry=Polygon.from_bounds(*info.bounds),
+                        properties=info,
                     )
 
     ############################################################################