Merge pull request #214 from csiro-coasts/centroid_coordinates

Add Grid.centroid_coordinates attribute

Author: mx-moth

docs/releases/development.rst

@@ -11,3 +11,5 @@ Next release (in development)
 * Add example showing
   :ref:`multiple ways of plotting vector data <sphx_glr_examples_plot-vector-methods.py>`
   (:pr:`213`).
+* Add :attr:`.Grid.centroid_coordinates` attribute
+  (:pr:`214`).

src/emsarray/conventions/_base.py

@@ -146,9 +146,36 @@ def centroid(self) -> numpy.ndarray:
         The centres of the geometry of this grid as a :class:`numpy.ndarray` of Shapely points.
         Defaults to the :func:`shapely.centroid` of :attr:`Grid.geometry`,
         but some conventions might have more specific ways of finding the centres.
+        If geometry is empty the corresponding index in this array will be None.
+
+        See also
+        ========
+        Grid.centroid
+            The same data represented as a two-dimensional array of coordinates.
+        Grid.mask
+            A boolean array indicating where the grid has geometry
         """
         return self.convention.make_geometry_centroid(self.grid_kind)
 
+    @cached_property
+    def centroid_coordinates(self) -> numpy.ndarray:
+        """
+        The centres of the geometry of this grid as a 2-dimensional :class:`numpy.ndarray`
+        with shape (:attr:`Grid.size`, 2).
+        If geometry is empty the corresponding row in this array will be `[numpy.nan, numpy.nan]`.
+
+        See also
+        ========
+        Grid.centroid
+            The same data represented as an array of :class:`shapely.Point`.
+        Grid.mask
+            A boolean array indicating where the grid has geometry
+        """
+        coordinates = numpy.full(fill_value=numpy.nan, shape=(self.size, 2))
+        _coordinates, indexes = shapely.get_coordinates(self.centroid, return_index=True)
+        coordinates[indexes] = _coordinates
+        return coordinates
+
     @abc.abstractmethod
     def ravel_index(self, index: Index) -> int:
         """
@@ -1492,15 +1519,12 @@ def polygons(self) -> numpy.ndarray:
     @cached_property
     @utils.deprecated(
         "dataset.ems.face_centres is deprecated. "
-        "Use dataset.ems.get_grid(data_array).centroid instead. "
+        "Use dataset.ems.get_grid(data_array).centroid_coordinates instead. "
         "For a list of coordinate pairs use shapely.get_coordinates(grid.centroid)."
     )
     def face_centres(self) -> numpy.ndarray:
         grid = self.grids[self.default_grid_kind]
-        centroid = grid.centroid
-        coords = numpy.full(fill_value=numpy.nan, shape=(grid.size, 2))
-        coords[centroid != None] = shapely.get_coordinates(centroid)  # noqa: E711
-        return cast(numpy.ndarray, coords)
+        return grid.centroid_coordinates
 
     @cached_property
     @utils.deprecated(

src/emsarray/plot/artists.py

@@ -255,10 +255,8 @@ def from_grid(
     def make_triangulation(grid: 'conventions.Grid') -> Triangulation:
         convention = grid.convention
         # Compute the Delaunay triangulation of the face centres
-        face_centres = grid.centroid
-        coords = numpy.full(fill_value=numpy.nan, shape=(grid.size, 2))
-        coords[face_centres != None] = shapely.get_coordinates(face_centres)  # noqa: E711
-        triangulation = Triangulation(coords[:, 0], coords[:, 1])
+        centres = grid.centroid_coordinates
+        triangulation = Triangulation(centres[:, 0], centres[:, 1])
 
         # Mask out any Triangles that are not contained within the dataset geometry.
         # These are either in concave areas of the geometry (e.g. an inlet or bay)
@@ -393,8 +391,7 @@ def from_grid(
         if not issubclass(grid.geometry_type, shapely.Polygon):
             raise ValueError("Grid must have polygon geometry")
 
-        coords = numpy.full(fill_value=numpy.nan, shape=(grid.size, 2))
-        coords[grid.centroid != None] = shapely.get_coordinates(grid.centroid)  # noqa: E711
+        coords = grid.centroid_coordinates
 
         # A Quiver needs some values when being initialized.
         # We don't always want to provide values to the quiver,

tests/conventions/test_base.py

@@ -536,3 +536,23 @@ def test_centroid():
     assert len(face_centres) == len(polygons)
     assert polygons[i] is None
     assert face_centres[i] is None
+
+
+def test_centroid_coordinates():
+    y_size, x_size = 10, 20
+    dataset = xarray.Dataset({
+        'temp': (['t', 'z', 'y', 'x'], numpy.random.standard_normal((5, 5, y_size, x_size))),
+        'botz': (['y', 'x'], numpy.random.standard_normal((y_size, x_size)) - 10),
+    })
+    convention = SimpleConvention(dataset)
+
+    grid = convention.grids['face']
+    xx, yy = numpy.meshgrid(
+        numpy.arange(x_size) + 0.5,
+        numpy.arange(y_size) + 0.5,
+    )
+    expected = numpy.c_[xx.flatten(), yy.flatten()]
+    expected[~grid.mask] = [numpy.nan, numpy.nan]
+    actual = grid.centroid_coordinates
+
+    numpy.testing.assert_equal(expected, actual)

tests/conventions/test_ugrid.py

@@ -463,34 +463,32 @@ def test_node_grid() -> None:
         assert node.equals_exact(shapely.Point([node_x, node_y]), 1e-6)
 
 
-def test_face_centres_from_variables():
+def test_face_centres_from_variables() -> None:
     dataset = make_dataset(width=3, make_face_coordinates=True)
     convention: UGrid = dataset.ems
-    face_grid = dataset.ems.grids['face']
+    face_grid = convention.grids['face']
 
-    face_centres = convention.default_grid.centroid
+    coordinates = face_grid.centroid_coordinates
     lons = dataset['Mesh2_face_x'].values
     lats = dataset['Mesh2_face_y'].values
     for face in range(dataset.sizes['nMesh2_face']):
         lon = lons[face]
         lat = lats[face]
         linear_index = face_grid.ravel_index((UGridKind.face, face))
-        point = face_centres[linear_index]
-        numpy.testing.assert_equal([point.x, point.y], [lon, lat])
+        numpy.testing.assert_equal(coordinates[linear_index], [lon, lat])
 
 
-def test_face_centres_from_centroids():
+def test_face_centres_from_centroids() -> None:
     dataset = make_dataset(width=3, make_face_coordinates=False)
     convention: UGrid = dataset.ems
-    face_grid = dataset.ems.grids['face']
-    face_centres = face_grid.centroid
+    face_grid = convention.grids['face']
+    coordinates = face_grid.centroid_coordinates
 
     for face in range(dataset.sizes['nMesh2_face']):
         linear_index = convention.ravel_index((UGridKind.face, face))
         polygon = face_grid.geometry[linear_index]
         lon, lat = polygon.centroid.coords[0]
-        point = face_centres[linear_index]
-        numpy.testing.assert_equal([point.x, point.y], [lon, lat])
+        numpy.testing.assert_equal(coordinates[linear_index], [lon, lat])
 
 
 def test_bounds(datasets: pathlib.Path):

tests/utils.py

@@ -456,8 +456,8 @@ def plot_geometry(
 
     dataset.ems.plot_geometry(axes)
     grid = dataset.ems.default_grid
-    centroid = shapely.get_coordinates(grid.centroid)
-    axes.scatter(centroid[:, 0], centroid[:, 1], c='red')
+    x, y = grid.centroid_coordinates.T
+    axes.scatter(x, y, c='red')
 
     if title is not None:
         axes.set_title(title)