Feat/algorithm in statistics (#726)

* enable algorithm in statistics endpoints

* update docs

CHANGES.md

@@ -1,5 +1,13 @@
 # Release Notes
 
+## 0.15.5 (2023-11-09)
+
+### titiler.core
+
+* add `algorithm` options for `/statistics` endpoints
+
+* switch from `BaseReader.statistics()` method to a combination of `BaseReader.preview()` and `ImageData.statistics()` methods to get the statistics
+
 ## 0.15.4 (2023-11-06)
 
 ### titiler.core

docs/src/endpoints/cog.md

@@ -54,7 +54,8 @@ The `/cog` routes are based on `titiler.core.factory.TilerFactory` but with `cog
     - **colormap_name** (str): rio-tiler color map name.
     - **return_mask** (bool): Add mask to the output data. Default is True.
     - **buffer** (float): Buffer on each side of the given tile. It must be a multiple of `0.5`. Output **tilesize** will be expanded to `tilesize + 2 * buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
-    - **padding** (int): Padding to apply to each tile edge. Helps reduce resampling artefacts along edges. Defaults to `0`.    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **padding** (int): Padding to apply to each tile edge. Helps reduce resampling artefacts along edges. Defaults to `0`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
     - **algorithm_params** (str): JSON encoded algorithm parameters.
 
 Example:
@@ -226,7 +227,8 @@ Example:
     - **colormap_name** (str): rio-tiler color map name.
     - **return_mask** (bool): Add mask to the output data. Default is True.
     - **buffer** (float): Buffer on each side of the given tile. It must be a multiple of `0.5`. Output **tilesize** will be expanded to `tilesize + 2 * buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
-    - **padding** (int): Padding to apply to each tile edge. Helps reduce resampling artefacts along edges. Defaults to `0`.    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **padding** (int): Padding to apply to each tile edge. Helps reduce resampling artefacts along edges. Defaults to `0`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
     - **algorithm_params** (str): JSON encoded algorithm parameters.
 
 Example:
@@ -261,7 +263,8 @@ Example:
     - **colormap_name** (str): rio-tiler color map name.
     - **return_mask** (bool): Add mask to the output data. Default is True.
     - **buffer** (float): Buffer on each side of the given tile. It must be a multiple of `0.5`. Output **tilesize** will be expanded to `tilesize + 2 * buffer` (e.g 0.5 = 257x257, 1.0 = 258x258).
-    - **padding** (int): Padding to apply to each tile edge. Helps reduce resampling artefacts along edges. Defaults to `0`.    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **padding** (int): Padding to apply to each tile edge. Helps reduce resampling artefacts along edges. Defaults to `0`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
     - **algorithm_params** (str): JSON encoded algorithm parameters.
 
 Example:
@@ -313,6 +316,8 @@ Advanced raster statistics
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
     - **resampling** (str): RasterIO resampling algorithm. Defaults to `nearest`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **algorithm_params** (str): JSON encoded algorithm parameters.
     - **categorical** (bool): Return statistics for categorical dataset, default is false.
     - **c** (array[float]): Pixels values for categories.
     - **p** (array[int]): Percentile values.
@@ -341,6 +346,8 @@ Example:
     - **unscale** (bool): Apply dataset internal Scale/Offset.
     - **resampling** (str): RasterIO resampling algorithm. Defaults to `nearest`.
     - **reproject** (str): WarpKernel resampling algorithm (only used when doing re-projection). Defaults to `nearest`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **algorithm_params** (str): JSON encoded algorithm parameters.
     - **categorical** (bool): Return statistics for categorical dataset, default is false.
     - **c** (array[float]): Pixels values for categories.
     - **p** (array[int]): Percentile values.

docs/src/endpoints/stac.md

@@ -382,6 +382,8 @@ Example:
     - **nodata** (str, int, float): Overwrite internal Nodata value.
     - **unscale** (bool): Apply dataset internal Scale/Offset.
     - **resampling** (str): RasterIO resampling algorithm. Defaults to `nearest`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **algorithm_params** (str): JSON encoded algorithm parameters.
     - **categorical** (bool): Return statistics for categorical dataset, default is false.
     - **c** (array[float]): Pixels values for categories.
     - **p** (array[int]): Percentile values.
@@ -413,6 +415,8 @@ Example:
     - **unscale** (bool): Apply dataset internal Scale/Offset.
     - **resampling** (str): RasterIO resampling algorithm. Defaults to `nearest`.
     - **reproject** (str): WarpKernel resampling algorithm (only used when doing re-projection). Defaults to `nearest`.
+    - **algorithm** (str): Custom algorithm name (e.g `hillshade`).
+    - **algorithm_params** (str): JSON encoded algorithm parameters.
     - **categorical** (bool): Return statistics for categorical dataset, default is false.
     - **c** (array[float]): Pixels values for categories.
     - **p** (array[int]): Percentile values.

src/titiler/core/tests/test_factories.py

@@ -525,6 +525,16 @@ def test_TilerFactory():
     assert min(resp["b1"]["histogram"][1]) == 5.0
     assert max(resp["b1"]["histogram"][1]) == 10.0
 
+    # Stats with Algorithm
+    response = client.get(
+        f"/statistics?url={DATA_DIR}/cog.tif&bidx=1&bidx=1&algorithm=normalizedIndex"
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/json"
+    resp = response.json()
+    assert len(resp) == 1
+    assert "(b1 - b1) / (b1 + b1)" in resp
+
     # POST - statistics
     response = client.post(
         f"/statistics?url={DATA_DIR}/cog.tif&bidx=1&bidx=1&bidx=1", json=feature
@@ -650,6 +660,18 @@ def test_TilerFactory():
     assert len(resp["properties"]["statistics"]["b1"]["histogram"][0]) == 4
     assert resp["properties"]["statistics"]["b1"]["histogram"][0][3] == 0
 
+    # Stats with Algorithm
+    response = client.post(
+        f"/statistics?url={DATA_DIR}/cog.tif&bidx=1&bidx=1&algorithm=normalizedIndex",
+        json=feature,
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/geo+json"
+    resp = response.json()
+    assert resp["type"] == "Feature"
+    assert len(resp["properties"]["statistics"]) == 1
+    assert "(b1 - b1) / (b1 + b1)" in resp["properties"]["statistics"]
+
     # Test with Algorithm
     response = client.get(f"/preview.tif?url={DATA_DIR}/dem.tif&return_mask=False")
     assert response.status_code == 200
@@ -865,6 +887,15 @@ def test_MultiBaseTilerFactory(rio):
     assert resp["B01_b1"]
     assert resp["B09_b1"]
 
+    # with Algorithm
+    response = client.get(
+        f"/statistics?url={DATA_DIR}/item.json&assets=B01&assets=B09&algorithm=normalizedIndex&asset_as_band=True"
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/json"
+    resp = response.json()
+    assert "(B09 - B01) / (B09 + B01)" in resp
+
     stac_feature = {
         "type": "FeatureCollection",
         "features": [
@@ -974,6 +1005,18 @@ def test_MultiBaseTilerFactory(rio):
     }
     assert props["B09_b1"]
 
+    # with Algorithm
+    response = client.post(
+        f"/statistics?url={DATA_DIR}/item.json&assets=B01&assets=B09&algorithm=normalizedIndex&asset_as_band=True",
+        json=stac_feature["features"][0],
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/geo+json"
+    resp = response.json()
+    props = resp["properties"]["statistics"]
+    assert len(props) == 1
+    assert "(B09 - B01) / (B09 + B01)" in props
+
 
 @attr.s
 class BandFileReader(MultiBandReader):
@@ -1142,6 +1185,15 @@ def test_MultiBandTilerFactory():
         "percentile_98",
     }
 
+    response = client.get(
+        f"/statistics?directory={DATA_DIR}&bands=B01&bands=B09&algorithm=normalizedIndex"
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/json"
+    resp = response.json()
+    assert len(resp) == 1
+    assert "(B09 - B01) / (B09 + B01)" in resp
+
     # POST - statistics
     band_feature = {
         "type": "FeatureCollection",
@@ -1277,6 +1329,17 @@ def test_MultiBandTilerFactory():
         "percentile_98",
     }
 
+    response = client.post(
+        f"/statistics?directory={DATA_DIR}&bands=B01&bands=B09&algorithm=normalizedIndex",
+        json=band_feature,
+    )
+    assert response.status_code == 200
+    assert response.headers["content-type"] == "application/geo+json"
+    resp = response.json()
+    props = resp["features"][0]["properties"]["statistics"]
+    assert len(props) == 1
+    assert "(B09 - B01) / (B09 + B01)" in props
+
     # default bands
     response = client.post(f"/statistics?directory={DATA_DIR}", json=band_feature)
     assert response.status_code == 200

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

@@ -407,6 +407,7 @@ def statistics(
             layer_params=Depends(self.layer_dependency),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_preview_dependency),
+            post_process=Depends(self.process_dependency),
             stats_params=Depends(self.stats_dependency),
             histogram_params=Depends(self.histogram_dependency),
             reader_params=Depends(self.reader_dependency),
@@ -415,10 +416,16 @@ def statistics(
             """Get Dataset statistics."""
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
-                    return src_dst.statistics(
+                    image = src_dst.preview(
                         **layer_params,
                         **image_params,
                         **dataset_params,
+                    )
+
+                    if post_process:
+                        image = post_process(image)
+
+                    return image.statistics(
                         **stats_params,
                         hist_options={**histogram_params},
                     )
@@ -447,6 +454,7 @@ def geojson_statistics(
             layer_params=Depends(self.layer_dependency),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_part_dependency),
+            post_process=Depends(self.process_dependency),
             stats_params=Depends(self.stats_dependency),
             histogram_params=Depends(self.histogram_dependency),
             reader_params=Depends(self.reader_dependency),
@@ -461,7 +469,7 @@ def geojson_statistics(
                 with self.reader(src_path, **reader_params) as src_dst:
                     for feature in fc:
                         shape = feature.model_dump(exclude_none=True)
-                        data = src_dst.feature(
+                        image = src_dst.feature(
                             shape,
                             shape_crs=coord_crs or WGS84_CRS,
                             dst_crs=dst_crs,
@@ -471,12 +479,15 @@ def geojson_statistics(
                         )
 
                         # Get the coverage % array
-                        coverage_array = data.get_coverage_array(
+                        coverage_array = image.get_coverage_array(
                             shape,
                             shape_crs=coord_crs or WGS84_CRS,
                         )
 
-                        stats = data.statistics(
+                        if post_process:
+                            image = post_process(image)
+
+                        stats = image.statistics(
                             **stats_params,
                             hist_options={**histogram_params},
                             coverage=coverage_array,
@@ -1225,6 +1236,7 @@ def statistics(
             layer_params=Depends(AssetsBidxExprParamsOptional),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_preview_dependency),
+            post_process=Depends(self.process_dependency),
             stats_params=Depends(self.stats_dependency),
             histogram_params=Depends(self.histogram_dependency),
             reader_params=Depends(self.reader_dependency),
@@ -1237,10 +1249,16 @@ def statistics(
                     if not layer_params.assets and not layer_params.expression:
                         layer_params.assets = src_dst.assets
 
-                    return src_dst.merged_statistics(
+                    image = src_dst.preview(
                         **layer_params,
                         **image_params,
                         **dataset_params,
+                    )
+
+                    if post_process:
+                        image = post_process(image)
+
+                    return image.statistics(
                         **stats_params,
                         hist_options={**histogram_params},
                     )
@@ -1268,6 +1286,7 @@ def geojson_statistics(
             dst_crs=Depends(DstCRSParams),
             layer_params=Depends(AssetsBidxExprParamsOptional),
             dataset_params=Depends(self.dataset_dependency),
+            post_process=Depends(self.process_dependency),
             image_params=Depends(self.img_part_dependency),
             stats_params=Depends(self.stats_dependency),
             histogram_params=Depends(self.histogram_dependency),
@@ -1286,7 +1305,7 @@ def geojson_statistics(
                         layer_params.assets = src_dst.assets
 
                     for feature in fc:
-                        data = src_dst.feature(
+                        image = src_dst.feature(
                             feature.model_dump(exclude_none=True),
                             shape_crs=coord_crs or WGS84_CRS,
                             dst_crs=dst_crs,
@@ -1295,7 +1314,10 @@ def geojson_statistics(
                             **dataset_params,
                         )
 
-                        stats = data.statistics(
+                        if post_process:
+                            image = post_process(image)
+
+                        stats = image.statistics(
                             **stats_params, hist_options={**histogram_params}
                         )
 
@@ -1416,6 +1438,7 @@ def statistics(
             bands_params=Depends(BandsExprParamsOptional),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_preview_dependency),
+            post_process=Depends(self.process_dependency),
             stats_params=Depends(self.stats_dependency),
             histogram_params=Depends(self.histogram_dependency),
             reader_params=Depends(self.reader_dependency),
@@ -1424,12 +1447,21 @@ def statistics(
             """Get Dataset statistics."""
             with rasterio.Env(**env):
                 with self.reader(src_path, **reader_params) as src_dst:
-                    return src_dst.statistics(
+                    # Default to all available bands
+                    if not bands_params.bands and not bands_params.expression:
+                        bands_params.bands = src_dst.bands
+
+                    image = src_dst.preview(
                         **bands_params,
                         **image_params,
                         **dataset_params,
-                        **stats_params,
-                        hist_options={**histogram_params},
+                    )
+
+                    if post_process:
+                        image = post_process(image)
+
+                    return image.statistics(
+                        **stats_params, hist_options={**histogram_params}
                     )
 
         # POST endpoint
@@ -1456,6 +1488,7 @@ def geojson_statistics(
             bands_params=Depends(BandsExprParamsOptional),
             dataset_params=Depends(self.dataset_dependency),
             image_params=Depends(self.img_part_dependency),
+            post_process=Depends(self.process_dependency),
             stats_params=Depends(self.stats_dependency),
             histogram_params=Depends(self.histogram_dependency),
             reader_params=Depends(self.reader_dependency),
@@ -1473,15 +1506,19 @@ def geojson_statistics(
                         bands_params.bands = src_dst.bands
 
                     for feature in fc:
-                        data = src_dst.feature(
+                        image = src_dst.feature(
                             feature.model_dump(exclude_none=True),
                             shape_crs=coord_crs or WGS84_CRS,
                             dst_crs=dst_crs,
                             **bands_params,
                             **image_params,
                             **dataset_params,
                         )
-                        stats = data.statistics(
+
+                        if post_process:
+                            image = post_process(image)
+
+                        stats = image.statistics(
                             **stats_params, hist_options={**histogram_params}
                         )