diff --git a/docs/releases/development.rst b/docs/releases/development.rst index 726f706..6b2b9fd 100644 --- a/docs/releases/development.rst +++ b/docs/releases/development.rst @@ -11,3 +11,5 @@ Next release (in development) * Add example showing :ref:`multiple ways of plotting vector data ` (:pr:`213`). +* Add :attr:`.Grid.centroid_coordinates` attribute + (:pr:`214`). diff --git a/src/emsarray/conventions/_base.py b/src/emsarray/conventions/_base.py index 24f01fb..c54e7ab 100644 --- a/src/emsarray/conventions/_base.py +++ b/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( diff --git a/src/emsarray/plot/artists.py b/src/emsarray/plot/artists.py index 51f7a2a..e85d1c9 100644 --- a/src/emsarray/plot/artists.py +++ b/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, diff --git a/tests/conventions/test_base.py b/tests/conventions/test_base.py index d575879..f674ae3 100644 --- a/tests/conventions/test_base.py +++ b/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) diff --git a/tests/conventions/test_ugrid.py b/tests/conventions/test_ugrid.py index e0363df..7b415fe 100644 --- a/tests/conventions/test_ugrid.py +++ b/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): diff --git a/tests/utils.py b/tests/utils.py index f0c69b4..fa85ac6 100644 --- a/tests/utils.py +++ b/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)