From 9ce42709d9b7333badfc77fc460c13f8374b6266 Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 30 Mar 2026 14:14:41 +0100 Subject: [PATCH 1/2] DistanceTo plugin modified to allow calculation of distance to features in angular sectors. Code changes made using copilot and some outputs validated but in need of more rigorous testing. --- .../generate_distance_to_feature.py | 238 +++++++++++++++++- .../generate_ancillaries/test_DistanceTo.py | 58 +++++ 2 files changed, 284 insertions(+), 12 deletions(-) diff --git a/improver/generate_ancillaries/generate_distance_to_feature.py b/improver/generate_ancillaries/generate_distance_to_feature.py index 2be93c1779..b26ae12021 100644 --- a/improver/generate_ancillaries/generate_distance_to_feature.py +++ b/improver/generate_ancillaries/generate_distance_to_feature.py @@ -5,13 +5,18 @@ """A module for generating a distance to feature ancillary cube.""" import os -from typing import List, Optional, Tuple +from typing import List, Optional, Tuple, Union +import iris +import numpy as np import pyproj from geopandas import GeoDataFrame, GeoSeries, clip +from iris.coords import AuxCoord, DimCoord from iris.cube import Cube from numpy import array, min, round -from shapely.geometry import Point +from shapely.geometry import Point, Polygon +from shapely.geometry.base import BaseGeometry +from shapely.ops import unary_union from improver import BasePlugin @@ -42,6 +47,7 @@ def __init__( self, epsg_projection: int, new_name: Optional[str] = None, + angle_pairs: Optional[List[Tuple[float, float]]] = None, buffer: float = 30000, clip_geometry_flag: bool = False, parallel: bool = False, @@ -59,6 +65,11 @@ def __init__( Areas projection suitable for the UK. new_name: The name of the output cube. + angle_pairs: + Optional list of angular sector bounds in degrees clockwise from + true north. Each element should be (start_angle, end_angle). + If provided, the shortest distance to geometry in each sector is + returned as an additional leading dimension on the output cube. buffer: A buffer distance in m. If the geometry is clipped, this distance will be added onto the outermost site locations to define the domain to clip @@ -78,11 +89,50 @@ def __init__( """ self.epsg_projection = epsg_projection self.new_name = new_name + self.angle_pairs = self.validate_angle_pairs(angle_pairs) self.buffer = buffer self.clip_geometry_flag = clip_geometry_flag self.parallel = parallel self.n_parallel_jobs = n_parallel_jobs + @staticmethod + def validate_angle_pairs( + angle_pairs: Optional[List[Tuple[float, float]]], + ) -> Optional[List[Tuple[float, float]]]: + """Validate sector angle pairs. + + Args: + angle_pairs: + Optional list of 2-tuples of angles in degrees. + + Returns: + Validated angle pairs. + + Raises: + ValueError: + If any angle pair is malformed. + """ + + if angle_pairs is None: + return None + + validated_pairs = [] + for pair in angle_pairs: + if len(pair) != 2: + raise ValueError( + "Each item in angle_pairs must contain exactly two angles." + ) + start_angle, end_angle = pair + if start_angle < 0 or start_angle > 360 or end_angle < 0 or end_angle > 360: + raise ValueError("Angles in angle_pairs must be between 0 and 360.") + if start_angle == end_angle and not (start_angle == 0 and end_angle == 360): + raise ValueError( + "Identical start and end angles are invalid unless using [0, 360]." + ) + validated_pairs.append((float(start_angle), float(end_angle))) + + return validated_pairs + @staticmethod def get_clip_values(points: List[float], buffer: float) -> List[float]: """Get the coordinates to use when clipping the geometry. This is determined by @@ -200,7 +250,64 @@ def project_geometry( site_points = site_points.to_crs(self.epsg_projection) return site_points, geometry_reprojection - def distance_to(self, site_points: GeoSeries, geometry: GeoDataFrame) -> List[int]: + @staticmethod + def _create_sector_geometry( + point: Point, + angle_pair: Tuple[float, float], + max_radius: float, + n_arc_points: int = 180, + ) -> Optional[BaseGeometry]: + """Create a sector polygon for the given site point. + + Args: + point: + Site location in projected coordinates. + angle_pair: + (start_angle, end_angle) in degrees clockwise from true north. + max_radius: + Radius used to create the sector polygon in metres. + n_arc_points: + Number of arc points used to approximate curved boundaries. + + Returns: + A Polygon representing the sector, or None for full-circle sectors. + """ + + start_angle, end_angle = angle_pair + if start_angle == 0 and end_angle == 360: + return None + + start_angle = start_angle % 360 + end_angle = end_angle % 360 + + sector_parts = ( + [(start_angle, end_angle)] + if end_angle > start_angle + else [(start_angle, 360.0), (0.0, end_angle)] + ) + + polygons = [] + for part_start, part_end in sector_parts: + angles = np.linspace(part_start, part_end, n_arc_points) + arc_points = [ + ( + point.x + max_radius * np.sin(np.deg2rad(angle)), + point.y + max_radius * np.cos(np.deg2rad(angle)), + ) + for angle in angles + ] + polygons.append( + Polygon([(point.x, point.y), *arc_points, (point.x, point.y)]) + ) + + return unary_union(polygons) + + def distance_to( + self, + site_points: GeoSeries, + geometry: GeoDataFrame, + angle_pair: Optional[Tuple[float, float]] = None, + ) -> List[float]: """Calculate the distance from each site point to the nearest feature in the geometry. @@ -209,11 +316,27 @@ def distance_to(self, site_points: GeoSeries, geometry: GeoDataFrame) -> List[in A GeoSeries containing the site points in the target projection. geometry: A GeoDataFrame containing the geometry in the target projection. + angle_pair: + Optional angular sector bounds, clockwise from true north. Returns: A list of distances from each site point to the nearest feature in the geometry rounded to the nearest metre.""" - def _distance_to_nearest(point: Point, geometry: GeoDataFrame) -> float: + bounds = geometry.total_bounds + corners = np.array( + [ + [bounds[0], bounds[1]], + [bounds[0], bounds[3]], + [bounds[2], bounds[1]], + [bounds[2], bounds[3]], + ] + ) + + def _distance_to_nearest( + point: Point, + geometry: GeoDataFrame, + angle_pair: Optional[Tuple[float, float]] = None, + ) -> float: """Calculate the distance from a point to the nearest feature in the geometry. Args: @@ -221,28 +344,54 @@ def _distance_to_nearest(point: Point, geometry: GeoDataFrame) -> float: A shapely Point object representing the site location. geometry: A GeoDataFrame containing the geometry in the target projection. + angle_pair: + Optional angular sector bounds, clockwise from true north. Returns: The distance from the point to the nearest feature in the geometry rounded to the nearest metre. """ - return round(min(point.distance(geometry.geometry))) + + geometry_to_search = geometry.geometry + if angle_pair is not None: + max_radius = ( + np.max(np.hypot(corners[:, 0] - point.x, corners[:, 1] - point.y)) + + 1.0 + ) + sector_polygon = self._create_sector_geometry( + point, angle_pair, max_radius=max_radius + ) + + if sector_polygon is not None: + intersected = geometry_to_search.intersection(sector_polygon) + geometry_to_search = intersected[ + ~intersected.is_empty & intersected.notnull() + ] + if geometry_to_search.empty: + return np.nan + + return round(min(point.distance(geometry_to_search))) if self.parallel: from joblib import Parallel, delayed parallel = Parallel(n_jobs=self.n_parallel_jobs, prefer="threads") output_generator = parallel( - delayed(_distance_to_nearest)(point, geometry) for point in site_points + delayed(_distance_to_nearest)(point, geometry, angle_pair) + for point in site_points ) distance_results = list(output_generator) else: distance_results = [] for point in site_points: - distance_results.append(_distance_to_nearest(point, geometry)) + distance_results.append( + _distance_to_nearest(point, geometry, angle_pair) + ) return distance_results - def create_output_cube(self, site_cube: Cube, data: List[int]) -> Cube: + def create_output_cube( + self, site_cube: Cube, data: Union[List[float], List[List[float]]] + ) -> Cube: """Create an output cube that will have the same metadata as the input site cube except the units are changed to metres and, if requested, the name of the output cube will be changed. @@ -252,13 +401,72 @@ def create_output_cube(self, site_cube: Cube, data: List[int]) -> Cube: The input cube containing site locations that are defined by latitude and longitude coordinates. data: - A list of distances from each site point to the nearest feature in the - geometry. + Distances from each site point to the nearest feature in the geometry, + optionally grouped by sector. Returns: A new cube containing the distances with the same metadata as input site cube but with updated units and name.""" - output_cube = site_cube.copy(data=array(data)) + if self.angle_pairs is None: + output_cube = site_cube.copy(data=array(data)) + else: + sector_data = array(data) + n_sectors = len(self.angle_pairs) + sector_coord = DimCoord( + np.arange(n_sectors, dtype=np.int32), long_name="sector", units="1" + ) + dim_coords_and_dims = [(sector_coord, 0)] + for coord in site_cube.dim_coords: + (dim,) = site_cube.coord_dims(coord) + dim_coords_and_dims.append((coord.copy(), dim + 1)) + + aux_coords_and_dims = [] + for coord in site_cube.aux_coords: + dims = site_cube.coord_dims(coord) + if dims: + aux_coords_and_dims.append( + (coord.copy(), tuple(dim + 1 for dim in dims)) + ) + else: + aux_coords_and_dims.append((coord.copy(), None)) + + angle_means = [] + for pair in self.angle_pairs: + if pair[0] > pair[1]: + angle_means.append((((pair[0] + pair[1]) / 2) - 180) % 360) + else: + angle_means.append(np.mean(pair)) + + aux_coords_and_dims.append( + ( + AuxCoord( + np.array( + angle_means, dtype=np.float32 + ), + bounds=np.array(self.angle_pairs, dtype=np.float32), + long_name="sector_angle_from_true_north", + units="degrees", + ), + 0, + ) + ) + + output_kwargs = { + "units": site_cube.units, + "attributes": site_cube.attributes.copy(), + "cell_methods": site_cube.cell_methods, + "dim_coords_and_dims": dim_coords_and_dims, + "aux_coords_and_dims": aux_coords_and_dims, + } + if site_cube.standard_name is not None: + output_kwargs["standard_name"] = site_cube.standard_name + else: + output_kwargs["long_name"] = site_cube.long_name + if site_cube.var_name is not None: + output_kwargs["var_name"] = site_cube.var_name + + output_cube = Cube(sector_data, **output_kwargs) + if self.new_name: output_cube.rename(self.new_name) output_cube.units = "m" @@ -316,6 +524,12 @@ def process(self, site_cube: Cube, geometry: GeoDataFrame) -> Cube: clipped_geometry = geometry_projection # Calculate the distance to the nearest feature in the geometry - distance_to_results = self.distance_to(site_coords, clipped_geometry) + if self.angle_pairs is None: + distance_to_results = self.distance_to(site_coords, clipped_geometry) + else: + distance_to_results = [ + self.distance_to(site_coords, clipped_geometry, angle_pair=angle_pair) + for angle_pair in self.angle_pairs + ] return self.create_output_cube(site_cube, distance_to_results) diff --git a/improver_tests/generate_ancillaries/test_DistanceTo.py b/improver_tests/generate_ancillaries/test_DistanceTo.py index 89014f29ec..b9003b6cdb 100644 --- a/improver_tests/generate_ancillaries/test_DistanceTo.py +++ b/improver_tests/generate_ancillaries/test_DistanceTo.py @@ -516,3 +516,61 @@ def test_distance_to_with_unsuitable_projection(single_site_cube, geometry_point ) with pytest.raises(ValueError, match=msg): DistanceTo(3112)(single_site_cube, geometry_point_laea) + + +def test_distance_to_with_angle_pairs(single_site_cube): + """Test that sector-based distances are returned with a leading sector dimension.""" + + single_site_cube.coord("latitude").points = 49.543481633 + single_site_cube.coord("longitude").points = -1.387510304 + + data = [ + Point(3500500, 3000800), # north, 300 m + Point(3500900, 3000500), # east, 400 m + Point(3500500, 3000400), # south, 100 m + Point(3500300, 3000500), # west, 200 m + ] + geometry = GeoDataFrame(geometry=data, crs="EPSG:3035") + + angle_pairs = [(315, 45), (45, 135), (135, 225), (225, 315)] + output_cube = DistanceTo(3035, angle_pairs=angle_pairs)(single_site_cube, geometry) + + expected = np.array([[300], [400], [100], [200]]) + np.testing.assert_allclose(output_cube.data, expected) + assert output_cube.shape == (4, 1) + np.testing.assert_allclose(output_cube.coord("sector").points, [0, 1, 2, 3]) + np.testing.assert_allclose( + output_cube.coord("sector_angle_from_true_north").points, + [0, 90, 180, 270], + ) + np.testing.assert_allclose( + output_cube.coord("sector_angle_from_true_north").bounds, + [ + [315, 45], [45, 135], [135, 225], [225, 315] + ], + ) + + +def test_distance_to_with_angle_pairs_no_feature_in_sector(single_site_cube): + """Test sectors with no geometry intersection return NaN.""" + + single_site_cube.coord("latitude").points = 49.543481633 + single_site_cube.coord("longitude").points = -1.387510304 + + geometry = GeoDataFrame(geometry=[Point(3500900, 3000500)], crs="EPSG:3035") + + angle_pairs = [(315, 45), (45, 135)] + output_cube = DistanceTo(3035, angle_pairs=angle_pairs)(single_site_cube, geometry) + + assert np.isnan(output_cube.data[0, 0]) + assert output_cube.data[1, 0] == 400 + + +def test_distance_to_invalid_angle_pairs(single_site_cube, geometry_point_laea): + """Test invalid sector definitions raise a ValueError.""" + + with pytest.raises( + ValueError, + match=r"Identical start and end angles are invalid unless using \[0, 360\].", + ): + DistanceTo(3035, angle_pairs=[(90, 90)])(single_site_cube, geometry_point_laea) From d202438f17c20df6f0bf684879533b51439d7afa Mon Sep 17 00:00:00 2001 From: "benjamin.ayliffe" Date: Mon, 30 Mar 2026 16:38:12 +0100 Subject: [PATCH 2/2] Style fixes. --- .../generate_ancillaries/generate_distance_to_feature.py | 5 +---- improver_tests/generate_ancillaries/test_DistanceTo.py | 4 +--- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/improver/generate_ancillaries/generate_distance_to_feature.py b/improver/generate_ancillaries/generate_distance_to_feature.py index b26ae12021..4140d70dae 100644 --- a/improver/generate_ancillaries/generate_distance_to_feature.py +++ b/improver/generate_ancillaries/generate_distance_to_feature.py @@ -7,7 +7,6 @@ import os from typing import List, Optional, Tuple, Union -import iris import numpy as np import pyproj from geopandas import GeoDataFrame, GeoSeries, clip @@ -440,9 +439,7 @@ def create_output_cube( aux_coords_and_dims.append( ( AuxCoord( - np.array( - angle_means, dtype=np.float32 - ), + np.array(angle_means, dtype=np.float32), bounds=np.array(self.angle_pairs, dtype=np.float32), long_name="sector_angle_from_true_north", units="degrees", diff --git a/improver_tests/generate_ancillaries/test_DistanceTo.py b/improver_tests/generate_ancillaries/test_DistanceTo.py index b9003b6cdb..96bb03289b 100644 --- a/improver_tests/generate_ancillaries/test_DistanceTo.py +++ b/improver_tests/generate_ancillaries/test_DistanceTo.py @@ -545,9 +545,7 @@ def test_distance_to_with_angle_pairs(single_site_cube): ) np.testing.assert_allclose( output_cube.coord("sector_angle_from_true_north").bounds, - [ - [315, 45], [45, 135], [135, 225], [225, 315] - ], + [[315, 45], [45, 135], [135, 225], [225, 315]], )