Fix/Voronoi issue

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
 

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):