Added named_arrays.geometry.point_in_polygon() to test if a given point is inside a polygon specified by its vertices.

README.md

@@ -2,7 +2,8 @@
 
 [![tests](https://github.com/sun-data/named-arrays/actions/workflows/tests.yml/badge.svg)](https://github.com/sun-data/named-arrays/actions/workflows/tests.yml)
 [![codecov](https://codecov.io/gh/sun-data/named-arrays/graph/badge.svg?token=1GhdcsgwO0)](https://codecov.io/gh/sun-data/named-arrays)
-[![Ruff](https://github.com/sun-data/named-arrays/actions/workflows/ruff.yml/badge.svg?branch=main)](https://github.com/sun-data/named-arrays/actions/workflows/ruff.yml)[![Documentation Status](https://readthedocs.org/projects/named-arrays/badge/?version=latest)](https://named-arrays.readthedocs.io/en/latest/?badge=latest)
+[![Ruff](https://github.com/sun-data/named-arrays/actions/workflows/ruff.yml/badge.svg?branch=main)](https://github.com/sun-data/named-arrays/actions/workflows/ruff.yml)
+[![Documentation Status](https://readthedocs.org/projects/named-arrays/badge/?version=latest)](https://named-arrays.readthedocs.io/en/latest/?badge=latest)
 [![PyPI version](https://badge.fury.io/py/named-arrays.svg)](https://badge.fury.io/py/named-arrays)
 
 `named-arrays` is an implementation of a [named tensor](https://nlp.seas.harvard.edu/NamedTensor), which assigns names to each axis of an n-dimensional array such as a numpy array.

named_arrays/__init__.py

@@ -35,6 +35,7 @@
 from . import ndfilters
 from . import colorsynth
 from . import numexpr
+from . import geometry
 from ._core import (
     QuantityLike,
     StartT,
@@ -402,6 +403,7 @@
     "ndfilters",
     "colorsynth",
     "numexpr",
+    "geometry",
     "QuantityLike",
     "StartT",
     "StopT",

named_arrays/_scalars/scalar_named_array_functions.py

@@ -14,6 +14,7 @@
 import regridding
 import named_arrays as na
 from . import scalars
+from ..geometry._point_in_polygon import _point_in_polygon_quantity
 
 __all__ = [
     "ASARRAY_LIKE_FUNCTIONS",
@@ -2082,3 +2083,53 @@ def evaluate(
         ndarray=result,
         axes=axes,
     )
+
+
+@_implements(na.geometry.point_in_polygon)
+def point_in_polygon(
+    x: na.AbstractScalarArray,
+    y: na.AbstractScalarArray,
+    vertices_x: na.AbstractScalarArray,
+    vertices_y: na.AbstractScalarArray,
+    axis: str,
+) -> na.ScalarArray:
+
+    try:
+        x = scalars._normalize(x)
+        y = scalars._normalize(y)
+        vertices_x = scalars._normalize(vertices_x)
+        vertices_y = scalars._normalize(vertices_y)
+    except scalars.ScalarTypeError:  # pragma: nocover
+        return NotImplemented
+
+    shape_vertices = na.shape_broadcasted(
+        x,
+        y,
+        vertices_x,
+        vertices_y,
+    )
+
+    shape = {a: shape_vertices[a] for a in shape_vertices if a != axis}
+
+    num_vertices = shape_vertices[axis]
+    shape_vertices = shape.copy()
+    shape_vertices[axis] = num_vertices
+
+    x = na.broadcast_to(x, shape)
+    y = na.broadcast_to(y, shape)
+    vertices_x = na.broadcast_to(vertices_x, shape_vertices)
+    vertices_y = na.broadcast_to(vertices_y, shape_vertices)
+
+    result = _point_in_polygon_quantity(
+        x=x.ndarray,
+        y=y.ndarray,
+        vertices_x=vertices_x.ndarray,
+        vertices_y=vertices_y.ndarray,
+    )
+
+    result = na.ScalarArray(
+        ndarray=result,
+        axes=tuple(shape),
+    )
+
+    return result

named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py

@@ -1343,3 +1343,42 @@ def evaluate(
             **kwargs,
         )
     )
+
+
+@_implements(na.geometry.point_in_polygon)
+def point_in_polygon(
+    x: na.AbstractScalar,
+    y: na.AbstractScalar,
+    vertices_x: na.AbstractScalar,
+    vertices_y: na.AbstractScalar,
+    axis: str,
+) -> na.UncertainScalarArray:
+
+    try:
+        x = uncertainties._normalize(x)
+        y = uncertainties._normalize(y)
+        vertices_x = uncertainties._normalize(vertices_x)
+        vertices_y = uncertainties._normalize(vertices_y)
+    except uncertainties.UncertainScalarTypeError:  # pragma: nocover
+        return NotImplemented
+
+    result_nominal = na.geometry.point_in_polygon(
+        x=x.nominal,
+        y=y.nominal,
+        vertices_x=vertices_x.nominal,
+        vertices_y=vertices_y.nominal,
+        axis=axis,
+    )
+
+    result_distribution = na.geometry.point_in_polygon(
+        x=x.distribution,
+        y=y.distribution,
+        vertices_x=vertices_x.distribution,
+        vertices_y=vertices_y.distribution,
+        axis=axis,
+    )
+
+    return na.UncertainScalarArray(
+        nominal=result_nominal,
+        distribution=result_distribution,
+    )

named_arrays/geometry/__init__.py

@@ -0,0 +1,7 @@
+"""Computational geometry routines."""
+
+from ._point_in_polygon import point_in_polygon
+
+__all__ = [
+    "point_in_polygon",
+]

named_arrays/geometry/_point_in_polygon.py

@@ -0,0 +1,191 @@
+from typing import TypeVar
+import numpy as np
+import astropy.units as u
+import numba
+import regridding
+import named_arrays as na
+
+PointT = TypeVar("PointT", bound="float | u.Quantity | na.AbstractScalar")
+VertexT = TypeVar("VertexT", bound="na.AbstractScalar")
+
+def point_in_polygon(
+    x: PointT,
+    y: PointT,
+    vertices_x: VertexT,
+    vertices_y: VertexT,
+    axis: str,
+) -> PointT | VertexT:
+    """
+    Check if a given point is inside or on the boundary of a polygon.
+
+    This function is a wrapper around
+    :func:`regridding.geometry.point_is_inside_polygon`.
+
+    Parameters
+    ----------
+    x
+        The :math:`x`-coordinates of the test points.
+    y
+        The :math:`y`-coordinates of the test points.
+    vertices_x
+        The :math:`x`-coordinates of the polygon's vertices.
+    vertices_y
+        The :math:`y`-coordinates of the polygon's vertices.
+    axis
+        The logical axis representing the different vertices of the polygon.
+
+    Examples
+    --------
+
+    Check if some random points are inside a randomly-generated polygon.
+
+    .. jupyter-execute::
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        import named_arrays as na
+
+        # Define a random polygon
+        axis = "vertex"
+        num_vertices = 7
+        radius = na.random.uniform(5, 15, shape_random={axis: num_vertices})
+        angle = na.linspace(0, 2 * np.pi, axis=axis, num=num_vertices)
+        vertices_x = radius * np.cos(angle)
+        vertices_y = radius * np.sin(angle)
+
+        # Define some random points
+        x = na.random.uniform(-20, 20, shape_random=dict(r=1000))
+        y = na.random.uniform(-20, 20, shape_random=dict(r=1000))
+
+        # Select which points are inside the polygon
+        where = na.geometry.point_in_polygon(
+            x=x,
+            y=y,
+            vertices_x=vertices_x,
+            vertices_y=vertices_y,
+            axis=axis,
+        )
+
+        # Plot the results as a scatter plot
+        fig, ax = plt.subplots()
+        na.plt.fill(
+            vertices_x,
+            vertices_y,
+            ax=ax,
+            facecolor="none",
+            edgecolor="black",
+        )
+        na.plt.scatter(
+            x,
+            y,
+            where=where,
+            ax=ax,
+        );
+    """
+    return na._named_array_function(
+        func=point_in_polygon,
+        x=x,
+        y=y,
+        vertices_x=vertices_x,
+        vertices_y=vertices_y,
+        axis=axis,
+    )
+
+
+def _point_in_polygon_quantity(
+    x: u.Quantity,
+    y: u.Quantity,
+    vertices_x: u.Quantity,
+    vertices_y: u.Quantity,
+) -> np.ndarray:
+    """
+    Check if a given point is inside or on the boundary of a polygon.
+
+    Parameters
+    ----------
+    x
+        The :math:`x`-coordinates of the test points.
+    y
+        The :math:`y`-coordinates of the test points.
+    vertices_x
+        The :math:`x`-coordinates of the polygon's vertices.
+        The last axis should represent the different vertices of the polygon.
+    vertices_y
+        The :math:`y`-coordinates of the polygon's vertices.
+        The last axis should represent the different vertices of the polygon.
+    """
+
+    if isinstance(x, u.Quantity):
+        unit = x.unit
+        y = y.to_value(unit)
+        vertices_x = vertices_x.to_value(unit)
+        vertices_y = vertices_y.to_value(unit)
+
+    shape_points = np.broadcast(x, y).shape
+    shape_vertices = np.broadcast(vertices_x, vertices_y).shape
+
+    num_vertices = shape_vertices[~0]
+
+    shape_points = np.broadcast_shapes(shape_points, shape_vertices[:~0])
+    shape_vertices = shape_points + (num_vertices,)
+
+    x = np.broadcast_to(x, shape_points)
+    y = np.broadcast_to(y, shape_points)
+
+    vertices_x = np.broadcast_to(vertices_x, shape_vertices)
+    vertices_y = np.broadcast_to(vertices_y, shape_vertices)
+
+    result = _point_in_polygon_numba(
+        x=x.reshape(-1),
+        y=y.reshape(-1),
+        vertices_x=vertices_x.reshape(-1, num_vertices),
+        vertices_y=vertices_y.reshape(-1, num_vertices),
+    )
+
+    result = result.reshape(shape_points)
+
+    return result
+
+@numba.njit(cache=True, parallel=True)
+def _point_in_polygon_numba(
+    x: np.ndarray,
+    y: np.ndarray,
+    vertices_x: np.ndarray,
+    vertices_y: np.ndarray,
+) -> np.ndarray:  # pragma: nocover
+    """
+    Numba-accelerated check if a given point is inside or on the boundary of a polygon.
+
+    Vectorized version of :func:`regridding.geometry.point_is_inside_polygon`.
+
+    Parameters
+    ----------
+    x
+        The :math:`x`-coordinates of the test points.
+        Should be 1-dimensional.
+    y
+        The :math:`y`-coordinates of the test points.
+        Should be 1-dimensional, with the same number of elements as `x`.
+    vertices_x
+        The :math:`x`-coordinates of the polygon's vertices.
+        Should be 2-dimensional, where the first axis has the same number
+        of elements as `x`.
+    vertices_y
+        The :math:`y`-coordinates of the polygon's vertices.
+        Should be 2-dimensional, where the last axis has the same number
+        of elements as `vertices_y`.
+    """
+
+    num_pts, num_vertices = vertices_x.shape
+
+    result = np.empty(num_pts, dtype=np.bool)
+
+    for i in numba.prange(num_pts):
+        result[i] = regridding.geometry.point_is_inside_polygon(
+            x=x[i],
+            y=y[i],
+            vertices_x=vertices_x[i],
+            vertices_y=vertices_y[i],
+        )
+
+    return result