From 6477185b1a33d690cd489b07bb5efc8a97c6e719 Mon Sep 17 00:00:00 2001 From: Joe Moorhouse <5102656+joemoorhouse@users.noreply.github.com> Date: Sun, 6 Oct 2024 20:19:16 +0100 Subject: [PATCH] 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 | 45 ++++++++++++++++++----------- src/physrisk/kernel/hazard_model.py | 16 +++++----- tests/data/events_retrieval_test.py | 21 ++++++++++++-- 3 files changed, 56 insertions(+), 26 deletions(-) diff --git a/src/physrisk/data/zarr_reader.py b/src/physrisk/data/zarr_reader.py index e7315e29..6d8a3f7f 100644 --- a/src/physrisk/data/zarr_reader.py +++ b/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,23 +299,26 @@ 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. @@ -316,13 +326,14 @@ def get_max_curves_on_grid( 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 diff --git a/src/physrisk/kernel/hazard_model.py b/src/physrisk/kernel/hazard_model.py index 235c7404..86e84f4f 100644 --- a/src/physrisk/kernel/hazard_model.py +++ b/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 diff --git a/tests/data/events_retrieval_test.py b/tests/data/events_retrieval_test.py index ef88adc3..ba7ab565 100644 --- a/tests/data/events_retrieval_test.py +++ b/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] + )