|
| 1 | +from typing import TypeVar |
| 2 | +import numpy as np |
| 3 | +import astropy.units as u |
| 4 | +import numba |
| 5 | +import regridding |
| 6 | +import named_arrays as na |
| 7 | + |
| 8 | +PointT = TypeVar("PointT", bound="float | u.Quantity | na.AbstractScalar") |
| 9 | +VertexT = TypeVar("VertexT", bound="na.AbstractScalar") |
| 10 | + |
| 11 | +def point_in_polygon( |
| 12 | + x: PointT, |
| 13 | + y: PointT, |
| 14 | + vertices_x: VertexT, |
| 15 | + vertices_y: VertexT, |
| 16 | + axis: str, |
| 17 | +) -> PointT | VertexT: |
| 18 | + """ |
| 19 | + Check if a given point is inside or on the boundary of a polygon. |
| 20 | +
|
| 21 | + This function is a wrapper around |
| 22 | + :func:`regridding.geometry.point_is_inside_polygon`. |
| 23 | +
|
| 24 | + Parameters |
| 25 | + ---------- |
| 26 | + x |
| 27 | + The :math:`x`-coordinates of the test points. |
| 28 | + y |
| 29 | + The :math:`y`-coordinates of the test points. |
| 30 | + vertices_x |
| 31 | + The :math:`x`-coordinates of the polygon's vertices. |
| 32 | + vertices_y |
| 33 | + The :math:`y`-coordinates of the polygon's vertices. |
| 34 | + axis |
| 35 | + The logical axis representing the different vertices of the polygon. |
| 36 | +
|
| 37 | + Examples |
| 38 | + -------- |
| 39 | +
|
| 40 | + Check if some random points are inside a randomly-generated polygon. |
| 41 | +
|
| 42 | + .. jupyter-execute:: |
| 43 | +
|
| 44 | + import numpy as np |
| 45 | + import matplotlib.pyplot as plt |
| 46 | + import named_arrays as na |
| 47 | +
|
| 48 | + # Define a random polygon |
| 49 | + axis = "vertex" |
| 50 | + num_vertices = 7 |
| 51 | + radius = na.random.uniform(5, 15, shape_random={axis: num_vertices}) |
| 52 | + angle = na.linspace(0, 2 * np.pi, axis=axis, num=num_vertices) |
| 53 | + vertices_x = radius * np.cos(angle) |
| 54 | + vertices_y = radius * np.sin(angle) |
| 55 | +
|
| 56 | + # Define some random points |
| 57 | + x = na.random.uniform(-20, 20, shape_random=dict(r=1000)) |
| 58 | + y = na.random.uniform(-20, 20, shape_random=dict(r=1000)) |
| 59 | +
|
| 60 | + # Select which points are inside the polygon |
| 61 | + where = na.geometry.point_in_polygon( |
| 62 | + x=x, |
| 63 | + y=y, |
| 64 | + vertices_x=vertices_x, |
| 65 | + vertices_y=vertices_y, |
| 66 | + axis=axis, |
| 67 | + ) |
| 68 | +
|
| 69 | + # Plot the results as a scatter plot |
| 70 | + fig, ax = plt.subplots() |
| 71 | + na.plt.fill( |
| 72 | + vertices_x, |
| 73 | + vertices_y, |
| 74 | + ax=ax, |
| 75 | + facecolor="none", |
| 76 | + edgecolor="black", |
| 77 | + ) |
| 78 | + na.plt.scatter( |
| 79 | + x, |
| 80 | + y, |
| 81 | + where=where, |
| 82 | + ax=ax, |
| 83 | + ); |
| 84 | + """ |
| 85 | + return na._named_array_function( |
| 86 | + func=point_in_polygon, |
| 87 | + x=x, |
| 88 | + y=y, |
| 89 | + vertices_x=vertices_x, |
| 90 | + vertices_y=vertices_y, |
| 91 | + axis=axis, |
| 92 | + ) |
| 93 | + |
| 94 | + |
| 95 | +def _point_in_polygon_quantity( |
| 96 | + x: u.Quantity, |
| 97 | + y: u.Quantity, |
| 98 | + vertices_x: u.Quantity, |
| 99 | + vertices_y: u.Quantity, |
| 100 | +) -> np.ndarray: |
| 101 | + """ |
| 102 | + Check if a given point is inside or on the boundary of a polygon. |
| 103 | +
|
| 104 | + Parameters |
| 105 | + ---------- |
| 106 | + x |
| 107 | + The :math:`x`-coordinates of the test points. |
| 108 | + y |
| 109 | + The :math:`y`-coordinates of the test points. |
| 110 | + vertices_x |
| 111 | + The :math:`x`-coordinates of the polygon's vertices. |
| 112 | + The last axis should represent the different vertices of the polygon. |
| 113 | + vertices_y |
| 114 | + The :math:`y`-coordinates of the polygon's vertices. |
| 115 | + The last axis should represent the different vertices of the polygon. |
| 116 | + """ |
| 117 | + |
| 118 | + if isinstance(x, u.Quantity): |
| 119 | + unit = x.unit |
| 120 | + y = y.to_value(unit) |
| 121 | + vertices_x = vertices_x.to_value(unit) |
| 122 | + vertices_y = vertices_y.to_value(unit) |
| 123 | + |
| 124 | + shape_points = np.broadcast(x, y).shape |
| 125 | + shape_vertices = np.broadcast(vertices_x, vertices_y).shape |
| 126 | + |
| 127 | + num_vertices = shape_vertices[~0] |
| 128 | + |
| 129 | + shape_points = np.broadcast_shapes(shape_points, shape_vertices[:~0]) |
| 130 | + shape_vertices = shape_points + (num_vertices,) |
| 131 | + |
| 132 | + x = np.broadcast_to(x, shape_points) |
| 133 | + y = np.broadcast_to(y, shape_points) |
| 134 | + |
| 135 | + vertices_x = np.broadcast_to(vertices_x, shape_vertices) |
| 136 | + vertices_y = np.broadcast_to(vertices_y, shape_vertices) |
| 137 | + |
| 138 | + result = _point_in_polygon_numba( |
| 139 | + x=x.reshape(-1), |
| 140 | + y=y.reshape(-1), |
| 141 | + vertices_x=vertices_x.reshape(-1, num_vertices), |
| 142 | + vertices_y=vertices_y.reshape(-1, num_vertices), |
| 143 | + ) |
| 144 | + |
| 145 | + result = result.reshape(shape_points) |
| 146 | + |
| 147 | + return result |
| 148 | + |
| 149 | +@numba.njit(cache=True, parallel=True) |
| 150 | +def _point_in_polygon_numba( |
| 151 | + x: np.ndarray, |
| 152 | + y: np.ndarray, |
| 153 | + vertices_x: np.ndarray, |
| 154 | + vertices_y: np.ndarray, |
| 155 | +) -> np.ndarray: # pragma: nocover |
| 156 | + """ |
| 157 | + Numba-accelerated check if a given point is inside or on the boundary of a polygon. |
| 158 | +
|
| 159 | + Vectorized version of :func:`regridding.geometry.point_is_inside_polygon`. |
| 160 | +
|
| 161 | + Parameters |
| 162 | + ---------- |
| 163 | + x |
| 164 | + The :math:`x`-coordinates of the test points. |
| 165 | + Should be 1-dimensional. |
| 166 | + y |
| 167 | + The :math:`y`-coordinates of the test points. |
| 168 | + Should be 1-dimensional, with the same number of elements as `x`. |
| 169 | + vertices_x |
| 170 | + The :math:`x`-coordinates of the polygon's vertices. |
| 171 | + Should be 2-dimensional, where the first axis has the same number |
| 172 | + of elements as `x`. |
| 173 | + vertices_y |
| 174 | + The :math:`y`-coordinates of the polygon's vertices. |
| 175 | + Should be 2-dimensional, where the last axis has the same number |
| 176 | + of elements as `vertices_y`. |
| 177 | + """ |
| 178 | + |
| 179 | + num_pts, num_vertices = vertices_x.shape |
| 180 | + |
| 181 | + result = np.empty(num_pts, dtype=np.bool) |
| 182 | + |
| 183 | + for i in numba.prange(num_pts): |
| 184 | + result[i] = regridding.geometry.point_is_inside_polygon( |
| 185 | + x=x[i], |
| 186 | + y=y[i], |
| 187 | + vertices_x=vertices_x[i], |
| 188 | + vertices_y=vertices_y[i], |
| 189 | + ) |
| 190 | + |
| 191 | + return result |
0 commit comments