Add missing reprojection for ZarrReader max in shape. (#362)

* Add missing reprojection for ZarrReaders max in shape.

Signed-off-by: Joe Moorhouse <5102656+joemoorhouse@users.noreply.github.com>

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Signed-off-by: Joe Moorhouse <5102656+joemoorhouse@users.noreply.github.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

src/physrisk/data/zarr_reader.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import s3fs
+import shapely.ops
 import zarr
 from affine import Affine
 from pyproj import Transformer
@@ -180,9 +181,15 @@ def get_max_curves(
 
         t = z.attrs["transform_mat3x3"]  # type: ignore
         transform = Affine(t[0], t[1], t[2], t[3], t[4], t[5])
+        crs = z.attrs.get("crs", "epsg:4326")
         units: str = z.attrs.get("units", "default")
 
+        if crs.lower() != "epsg:4236":
+            transproj = Transformer.from_crs("epsg:4326", crs, always_xy=True).transform
+            shapes = [shapely.ops.transform(transproj, shape) for shape in shapes]
+
         matrix = np.array(~transform).reshape(3, 3).transpose()[:, :-1].reshape(6)
+
         transformed_shapes = [
             affinity.affine_transform(shape, matrix) for shape in shapes
         ]
@@ -292,37 +299,41 @@ def get_max_curves_on_grid(
         n_grid=5,
     ):
         """Get maximal intensity curve for a grid around a given latitude and longitude coordinate pair.
-            It is almost equivalent to:
-                self.get_max_curves
-                (
-                    set_id,
-                    [
-                        Polygon(
-                            (
-                                (x - 0.5 * delta_deg, y - 0.5 * delta_deg),
-                                (x - 0.5 * delta_deg, y + 0.5 * delta_deg),
-                                (x + 0.5 * delta_deg, y + 0.5 * delta_deg),
-                                (x + 0.5 * delta_deg, y - 0.5 * delta_deg)
-                            )
-                        )
-                        for x, y in zip(longitudes, latitudes)
-                    ]
-                    interpolation
+        It is almost equivalent to:
+        ```
+        self.get_max_curves
+        (
+            set_id,
+            [
+                Polygon(
+                    (
+                        (x - 0.5 * delta_deg, y - 0.5 * delta_deg),
+                        (x - 0.5 * delta_deg, y + 0.5 * delta_deg),
+                        (x + 0.5 * delta_deg, y + 0.5 * delta_deg),
+                        (x + 0.5 * delta_deg, y - 0.5 * delta_deg)
+                    )
                 )
+                for x, y in zip(longitudes, latitudes)
+            ]
+            interpolation
+        )
+        ```
+
         Args:
             set_id: string or tuple representing data set, converted into path by path_provider.
             longitudes: list of longitudes.
             latitudes: list of latitudes.
             interpolation: interpolation method, "floor", "linear", "max" or "min".
             delta_km: linear distance in kilometres of the side of the square grid surrounding a given position.
             n_grid: number of grid points along the latitude and longitude dimensions used for
-            calculating the maximal value.
+                calculating the maximal value.
 
         Returns:
             curves_max: numpy array of maximum intensity on the grid for a given coordinate pair
             (no. coordinate pairs, no. return periods).
             return_periods: return periods in years.
         """
+
         kilometres_per_degree = 110.574
         delta_deg = delta_km / kilometres_per_degree
 

src/physrisk/kernel/hazard_model.py

@@ -30,13 +30,15 @@ def __init__(
         """Create HazardDataRequest.
 
         Args:
-            event_type: type of hazard.
-            longitude: required longitude.
-            latitude: required latitude.
-            model: model identifier.
-            scenario: identifier of scenario, e.g. rcp8p5 (RCP 8.5).
-            year: projection year, e.g. 2080.
-            buffer: delimitation of the area for the hazard data expressed in metres (within [0,1000]).
+            hazard_type (Type[Hazard]): Type of hazard.
+            longitude (float): Longitude.
+            latitude (float): Latitude.
+            indicator_id (str): Hazard indicator identifier, e.g. "flood_depth".
+            scenario (str): Scenario identifier, e.g. "SSP585", "RCP8p5".
+            year (int): Year for which data required.
+            hint (Optional[HazardDataHint], optional): Hint, typically providing the data set path.
+                Defaults to None.
+            buffer (Optional[int], optional): If not None applies a buffer around the point in metres, within [0, 1000m]. Defaults to None.
         """
         self.hazard_type = hazard_type
         self.longitude = longitude

tests/data/events_retrieval_test.py

@@ -267,8 +267,8 @@ def test_zarr_geomax(self):
         )
 
     def test_reproject(self):
-        """Test adding data in a non-ESPG-4326 coordinate reference system. Check that attribute
-        end in the correct convertion."""
+        """Test adding data in a non-ESPG-4326 coordinate reference system. Check that the round tip yields
+        the correct results."""
         mocker = ZarrStoreMocker()
         lons = [1.1, -0.31]
         lats = [47.0, 52.0]
@@ -302,3 +302,20 @@ def test_reproject(self):
         numpy.testing.assert_equal(
             next(iter(response.values())).intensities, [1.0, 2.0, 3.0]
         )
+        # run with a non-point shape also (hits different code path)
+        response2 = hazard_model.get_hazard_data(
+            [
+                HazardDataRequest(
+                    RiverineInundation,
+                    lons[0],
+                    lats[0],
+                    indicator_id="",
+                    scenario="",
+                    year=2050,
+                    buffer=10,
+                )
+            ]
+        )
+        numpy.testing.assert_equal(
+            next(iter(response2.values())).intensities, [1.0, 2.0, 3.0]
+        )