Merge pull request #1071 from xcube-dev/forman-1069-stats_for_non_wgs…

…84_crs

Make "/statistics" endpoint work for non-WGS84 grid mappings

CHANGES.md

@@ -4,6 +4,9 @@
 
 * The `time` query parameter of the `/statistics` endpoint of xcube server has 
    now been made optional. (#1066)
+* The `/statistics` endpoint now supports datasets using non-WGS84 grid systems, 
+  expanding its compatibility with a wider range of geospatial datasets.
+  (#1069)
 
 ## Changes in 1.7.0
 

test/webapi/statistics/test_routes.py

@@ -2,6 +2,8 @@
 # Permissions are hereby granted under the terms of the MIT License:
 # https://opensource.org/licenses/MIT.
 
+import json
+
 from ..helpers import RoutesTestCase
 
 
@@ -14,11 +16,32 @@ def get_config_filename(self) -> str:
 
     def test_fetch_post_statistics_ok(self):
         response = self.fetch(
-            "/statistics/demo/conc_chl?time=2017-01-16+10:09:21",
+            "/statistics/demo/conc_chl?time=2017-01-30+10:46:34",
             method="POST",
-            body='{"type": "Point", "coordinates": [1.768, 51.465]}',
+            body='{"type": "Point", "coordinates": [1.262, 50.243]}',
+        )
+        self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert parsed_data["result"]["count"] == 1
+        assert round(parsed_data["result"]["minimum"], 3) == 9.173
+        assert round(parsed_data["result"]["maximum"], 3) == 9.173
+        assert round(parsed_data["result"]["mean"], 3) == 9.173
+        assert round(parsed_data["result"]["deviation"], 3) == 0.0
+
+        response = self.fetch(
+            "/statistics/cog_local/band_1",
+            method="POST",
+            body='{"type": "Point", "coordinates": [-105.810, 35.771]}',
         )
         self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert parsed_data["result"]["count"] == 1
+        assert round(parsed_data["result"]["minimum"], 3) == 102.0
+        assert round(parsed_data["result"]["maximum"], 3) == 102.0
+        assert round(parsed_data["result"]["mean"], 3) == 102.0
+        assert round(parsed_data["result"]["deviation"], 3) == 0.0
 
     def test_fetch_post_statistics_missing_time_with_time_dimension_dataset(self):
         response = self.fetch(
@@ -30,15 +53,15 @@ def test_fetch_post_statistics_missing_time_with_time_dimension_dataset(self):
 
     def test_fetch_post_statistics_missing_time_without_time_dimension_dataset(self):
         response = self.fetch(
-            "/statistics/cog_local/band-1",
+            "/statistics/cog_local/band_1",
             method="POST",
             body='{"type": "Point", "coordinates": [-105.591, 35.751]}',
         )
         self.assertResponseOK(response)
 
     def test_fetch_post_statistics_with_time_without_time_dimension_dataset(self):
         response = self.fetch(
-            "/statistics/cog_local/band-1?time=2017-01-16+10:09:21",
+            "/statistics/cog_local/band_1?time=2017-01-16+10:09:21",
             method="POST",
             body='{"type": "Point", "coordinates": [-105.591, 35.751]}',
         )
@@ -52,50 +75,107 @@ def test_fetch_post_statistics_invalid_geometry(self):
         response = self.fetch(
             "/statistics/demo/conc_chl?time=2017-01-16+10:09:21",
             method="POST",
-            body="[1.768, 51.465]",
+            body="[1.768, 51.465, 11.542]",
+        )
+        self.assertBadRequestResponse(
+            response, "Invalid " "GeoJSON geometry encountered"
+        )
+        response = self.fetch(
+            "/statistics/demo/conc_chl?time=2017-01-16+10:09:21",
+            method="POST",
+            body="[1.768]",
         )
         self.assertBadRequestResponse(
             response, "Invalid " "GeoJSON geometry encountered"
         )
 
-    def test_fetch_get_statistics_ok(self):
+    def test_crs_conversion_post_statistics_with_coordinates_outside_bounds(self):
         response = self.fetch(
-            "/statistics/demo/conc_chl?"
-            "lat=1.786&lon=51.465&time=2017-01-16+10:09:21",
-            method="GET",
+            "/statistics/cog_local/band_1",
+            method="POST",
+            body='{"type": "Point", "coordinates": [-125.810, 35.771]}',
+        )
+        self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert parsed_data["result"]["count"] == 0
+
+    def test_crs_conversion_post_statistics_with_coordinates_inside_bounds(self):
+        response = self.fetch(
+            "/statistics/cog_local/band_1",
+            method="POST",
+            body='{"type": "Point", "coordinates": [-105.810, 35.171]}',
         )
         self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert parsed_data["result"]["count"] == 1
+        assert round(parsed_data["result"]["minimum"], 3) == 220.0
+        assert round(parsed_data["result"]["maximum"], 3) == 220.0
+        assert round(parsed_data["result"]["mean"], 3) == 220.0
+        assert round(parsed_data["result"]["deviation"], 3) == 0.0
 
     def test_fetch_get_statistics_missing_time_with_time_dimension_dataset(self):
         response = self.fetch(
-            "/statistics/demo/conc_chl?lat=1.786&lon=51.465", method="GET"
+            "/statistics/demo/conc_chl?lon=1.786&lat=51.465", method="GET"
         )
         self.assertBadRequestResponse(response, "Missing " "query parameter 'time'")
 
     def test_fetch_get_statistics_missing_time_without_time_dimension_dataset(self):
         response = self.fetch(
-            "/statistics/cog_local/band-1?lat=-105.591&" "lon=35.751&type=Point",
+            "/statistics/cog_local/band_1?lon=-105.591&" "lat=35.751&type=Point",
             method="GET",
         )
         self.assertResponseOK(response)
 
     def test_fetch_get_statistics_with_time_without_time_dimension_dataset(self):
         response = self.fetch(
-            "/statistics/cog_local/band-1?lat=-105.591&lon=35.751&"
+            "/statistics/cog_local/band_1?lon=-105.591&lat=35.751&"
             "type=Point&time=2017-01-16+10:09:21",
             method="GET",
-            body='{"type": "Point", "coordinates": [-105.591, 35.751]}',
         )
         self.assertBadRequestResponse(
             response,
             "Query parameter 'time' must not be given since "
             "dataset does not contain a 'time' dimension",
         )
 
-    def test_fetch_get_statistics_invalid_geometry(self):
+    def test_fetch_get_statistics(self):
+        response = self.fetch(
+            "/statistics/demo/conc_chl?time=2017-01-30+10:46:34&"
+            "lon=1.262&lat=50.243",
+            method="GET",
+        )
+        self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert round(parsed_data["result"]["value"], 3) == 9.173
+
+        response = self.fetch(
+            "/statistics/cog_local/band_1?lon=-105.810&lat=35.771&type=Point",
+            method="GET",
+        )
+        self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert round(parsed_data["result"]["value"], 3) == 102.0
+
+    def test_crs_conversion_get_statistics_with_coordinates_outside_bounds(self):
+        response = self.fetch(
+            "/statistics/cog_local/band_1?lon=-125.810&lat=35.771&type=Point",
+            method="GET",
+        )
+        self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert parsed_data["result"] == {}
+
+    def test_crs_conversion_get_statistics_with_coordinates_inside_bounds(self):
         response = self.fetch(
-            "/statistics/demo/conc_chl?time=2017-01-16+10:09:21&"
-            "lon=1.768&lat=51.465",
+            "/statistics/cog_local/band_1?lon=-105.810&lat=35.171&type=Point",
             method="GET",
         )
         self.assertResponseOK(response)
+        decoded_data = response.data.decode("utf-8")
+        parsed_data = json.loads(decoded_data)
+        assert round(parsed_data["result"]["value"], 3) == 220.0

xcube/webapi/statistics/controllers.py

@@ -1,12 +1,15 @@
-from collections.abc import Mapping
 from typing import Any, Union
 
 import numpy as np
-import xarray as xr
+import pyproj
 import shapely
+import shapely.ops
+import xarray as xr
 
+from xcube.constants import CRS_CRS84
 from xcube.constants import LOG
 from xcube.core.geom import get_dataset_geometry
+from xcube.core.geom import normalize_geometry
 from xcube.core.geom import mask_dataset_by_geometry
 from xcube.core.varexpr import VarExprContext
 from xcube.core.varexpr import split_var_assignment
@@ -70,9 +73,17 @@ def _compute_statistics(
     else:
         compact_mode = False
         try:
-            geometry = shapely.geometry.shape(geometry)
+            geometry = normalize_geometry(geometry)
         except (TypeError, ValueError, AttributeError) as e:
-            raise ApiError.BadRequest("Invalid GeoJSON geometry " "encountered") from e
+            raise ApiError.BadRequest(
+                f"Invalid GeoJSON geometry " "encountered: {e}"
+            ) from e
+
+    if geometry is not None and not grid_mapping.crs.is_geographic:
+        project = pyproj.Transformer.from_crs(
+            CRS_CRS84, grid_mapping.crs, always_xy=True
+        ).transform
+        geometry = shapely.ops.transform(project, geometry)
 
     nan_result = NAN_RESULT_COMPACT if compact_mode else NAN_RESULT