Merge pull request #156 from csiro-coasts/test-valid-polygons

Test polygons are valid after constructing them

Author: mx-moth

docs/releases/development.rst

@@ -6,3 +6,8 @@ Next release (in development)
   in CFGrid2D datasets with no cell bounds (:pr:`154`).
 * Improved speed of triangulation for convex polygons
   (:pr:`151`).
+* Check all polygons in a dataset are valid as part of generating them.
+  This will slow down opening new datasets slightly,
+  but the trade off is worth the added security
+  after the invalid polygons found in :pr:`154`
+  (:pr:`156`).

src/emsarray/conventions/_base.py

@@ -8,15 +8,15 @@
 from typing import TYPE_CHECKING, Any, Generic, Literal, TypeVar, cast
 
 import numpy
+import shapely
 import xarray
-from shapely import unary_union
 from shapely.geometry import MultiPolygon, Point, Polygon
 from shapely.geometry.base import BaseGeometry
 from shapely.strtree import STRtree
 
 from emsarray import utils
 from emsarray.compat.shapely import SpatialIndex
-from emsarray.exceptions import NoSuchCoordinateError
+from emsarray.exceptions import InvalidPolygonWarning, NoSuchCoordinateError
 from emsarray.operations import depth, point_extraction
 from emsarray.plot import (
     _requires_plot, animate_on_figure, make_plot_title, plot_on_figure,
@@ -1264,8 +1264,8 @@ def make_quiver(
 
         return Quiver(axes, x, y, *values, **kwargs)
 
-    @property
-    @abc.abstractmethod
+    @cached_property
+    @utils.timed_func
     def polygons(self) -> numpy.ndarray:
         """A :class:`numpy.ndarray` of :class:`shapely.Polygon` instances
         representing the cells in this dataset.
@@ -1286,6 +1286,24 @@ def polygons(self) -> numpy.ndarray:
         :meth:`ravel_index`
         :attr:`mask`
         """
+        polygons = self._make_polygons()
+
+        not_none = (polygons != None)  # noqa: E711
+        invalid_polygon_indices = numpy.flatnonzero(not_none & ~shapely.is_valid(polygons))
+        if len(invalid_polygon_indices):
+            indices_str = numpy.array2string(
+                invalid_polygon_indices, max_line_width=None, threshold=5)
+            warnings.warn(
+                f"Dropping invalid polygons at indices {indices_str}",
+                category=InvalidPolygonWarning)
+            polygons[invalid_polygon_indices] = None
+            not_none[invalid_polygon_indices] = False
+
+        polygons.flags.writeable = False
+        return polygons
+
+    @abc.abstractmethod
+    def _make_polygons(self) -> numpy.ndarray:
         pass
 
     @cached_property
@@ -1337,7 +1355,7 @@ def geometry(self) -> Polygon | MultiPolygon:
         This is equivalent to the union of all polygons in the dataset,
         although specific conventions may have a simpler way of constructing this.
         """
-        return unary_union(self.polygons[self.mask])
+        return shapely.unary_union(self.polygons[self.mask])
 
     @cached_property
     def bounds(self) -> Bounds:

src/emsarray/conventions/arakawa_c.py

@@ -260,9 +260,7 @@ def unpack_index(self, index: ArakawaCIndex) -> tuple[ArakawaCGridKind, Sequence
     def pack_index(self, grid_kind: ArakawaCGridKind, indexes: Sequence[int]) -> ArakawaCIndex:
         return cast(ArakawaCIndex, (grid_kind, *indexes))
 
-    @cached_property
-    @utils.timed_func
-    def polygons(self) -> numpy.ndarray:
+    def _make_polygons(self) -> numpy.ndarray:
         # Make an array of shape (j, i, 2) of all the nodes
         grid = numpy.stack([self.node.longitude.values, self.node.latitude.values], axis=-1)
 

src/emsarray/conventions/grid.py

@@ -408,9 +408,7 @@ def check_dataset(cls, dataset: xarray.Dataset) -> int | None:
 
         return Specificity.LOW
 
-    @cached_property
-    @utils.timed_func
-    def polygons(self) -> numpy.ndarray:
+    def _make_polygons(self) -> numpy.ndarray:
         lon_bounds = self.topology.longitude_bounds.values
         lat_bounds = self.topology.latitude_bounds.values
 
@@ -576,9 +574,7 @@ def check_dataset(cls, dataset: xarray.Dataset) -> int | None:
 
         return Specificity.LOW
 
-    @cached_property
-    @utils.timed_func
-    def polygons(self) -> numpy.ndarray:
+    def _make_polygons(self) -> numpy.ndarray:
         # Construct polygons from the bounds of the cells
         lon_bounds = self.topology.longitude_bounds.values
         lat_bounds = self.topology.latitude_bounds.values

src/emsarray/conventions/ugrid.py

@@ -1076,9 +1076,7 @@ def grid_kinds(self) -> frozenset[UGridKind]:
             items.append(UGridKind.edge)
         return frozenset(items)
 
-    @cached_property
-    @utils.timed_func
-    def polygons(self) -> numpy.ndarray:
+    def _make_polygons(self) -> numpy.ndarray:
         """Generate list of Polygons"""
         # X,Y coords of each node
         topology = self.topology

src/emsarray/exceptions.py

@@ -10,13 +10,19 @@ class EmsarrayError(Exception):
     """
 
 
+class EmsarrayWarning(Warning):
+    """
+    Base class for all emsarray-specific warning classes.
+    """
+
+
 class ConventionViolationError(EmsarrayError):
     """
     A dataset violates its conventions in a way that is not recoverable.
     """
 
 
-class ConventionViolationWarning(UserWarning):
+class ConventionViolationWarning(EmsarrayWarning):
     """
     A dataset violates its conventions in a way that we can handle.
     For example, an attribute has an invalid type,
@@ -30,3 +36,9 @@ class NoSuchCoordinateError(KeyError, EmsarrayError):
     such as in :attr:`.Convention.time_coordinate` and
     :attr:`.Convention.depth_coordinate`.
     """
+
+
+class InvalidPolygonWarning(EmsarrayWarning):
+    """
+    A polygon in a dataset was invalid or not simple.
+    """

src/emsarray/operations/triangulate.py

@@ -147,9 +147,6 @@ def _triangulate_polygon(polygon: Polygon) -> list[tuple[Vertex, Vertex, Vertex]
     :func:`triangulate_dataset`,
     `Polygon triangulation <https://en.wikipedia.org/wiki/Polygon_triangulation>`_
     """
-    if not polygon.is_simple:
-        raise ValueError("_triangulate_polygon only supports simple polygons")
-
     # The 'ear clipping' method used below is correct for all polygons, but not
     # performant. If the polygon is convex we can use a shortcut method.
     if polygon.equals(polygon.convex_hull):
@@ -174,9 +171,6 @@ def _triangulate_polygon(polygon: Polygon) -> list[tuple[Vertex, Vertex, Vertex]
     # Most polygons will be either squares, convex quadrilaterals, or convex
     # polygons.
 
-    # Maintain a consistent winding order
-    polygon = polygon.normalize()
-
     triangles: list[tuple[Vertex, Vertex, Vertex]] = []
     # Note that shapely polygons with n vertices will be closed, and thus have
     # n+1 coordinates. We trim that superfluous coordinate off in the next line

tests/conventions/test_base.py

@@ -118,8 +118,7 @@ def wind(
     def drop_geometry(self) -> xarray.Dataset:
         return self.dataset
 
-    @cached_property
-    def polygons(self) -> numpy.ndarray:
+    def _make_polygons(self) -> numpy.ndarray:
         height, width = self.shape
         # Each polygon is a box from (x, y, x+1, y+1),
         # however the polygons around the edge are masked out with None.