diff --git a/CHANGELOG.md b/CHANGELOG.md index 7743f30..16b9282 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,7 @@ Keep it human-readable, your future self will thank you! - Added `CutOutMask` class to create a mask for a cutout. (#68) - Added `MissingZarrVariable` and `NotMissingZarrVariable` classes to create a mask for missing zarr variables. (#68) - feat: Add CONTRIBUTORS.md file. (#72) +- Fixed issue when computing area weights with scipy.Voronoi. (#79) ### Changed diff --git a/src/anemoi/graphs/nodes/attributes.py b/src/anemoi/graphs/nodes/attributes.py index 1498e6d..4a7ad10 100644 --- a/src/anemoi/graphs/nodes/attributes.py +++ b/src/anemoi/graphs/nodes/attributes.py @@ -106,6 +106,8 @@ class AreaWeights(BaseNodeAttribute): Radius of the sphere. centre : np.ndarray Centre of the sphere. + fill_value : float + Value to fill the empty regions. Methods ------- @@ -118,11 +120,13 @@ def __init__( norm: str | None = None, radius: float = 1.0, centre: np.ndarray = np.array([0, 0, 0]), + fill_value: float = 0.0, dtype: str = "float32", ) -> None: super().__init__(norm, dtype) self.radius = radius self.centre = centre + self.fill_value = fill_value def get_raw_values(self, nodes: NodeStorage, **kwargs) -> np.ndarray: """Compute the area associated to each node. @@ -144,13 +148,19 @@ def get_raw_values(self, nodes: NodeStorage, **kwargs) -> np.ndarray: latitudes, longitudes = nodes.x[:, 0], nodes.x[:, 1] points = latlon_rad_to_cartesian((np.asarray(latitudes), np.asarray(longitudes))) sv = SphericalVoronoi(points, self.radius, self.centre) + mask = np.array([bool(i) for i in sv.regions]) + sv.regions = [region for region in sv.regions if region] + # compute the area weight without empty regions area_weights = sv.calculate_areas() + # add them back with zero weight + result = np.ones(points.shape[0]) * self.fill_value + result[mask] = area_weights LOGGER.debug( "There are %d of weights, which (unscaled) add up a total weight of %.2f.", - len(area_weights), - np.array(area_weights).sum(), + len(result), + np.array(result).sum(), ) - return area_weights + return result class BooleanBaseNodeAttribute(BaseNodeAttribute, ABC):