Coverage for python/lsst/images/_polygon.py: 91%
118 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-27 08:32 +0000
1# This file is part of lsst-images.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# Use of this source code is governed by a 3-clause BSD-style
10# license that can be found in the LICENSE file.
12from __future__ import annotations
14__all__ = ("Polygon", "Region", "RegionSerializationModel")
16from typing import TYPE_CHECKING, Any, Literal, overload
18import numpy as np
19import numpy.typing as npt
20import pydantic
21import shapely
23from ._geom import Bounds, Box
24from .utils import round_half_down, round_half_up
26if TYPE_CHECKING:
27 try:
28 from lsst.afw.geom import Polygon as LegacyPolygon
29 except ImportError:
30 type LegacyPolygon = Any # type: ignore[no-redef]
33class Region:
34 """A 2-d Euclidean region represented as one or more polygons with
35 optional holes.
37 Parameters
38 ----------
39 geometry
40 A polygon or multi-polygon from the Shapely library.
41 """
43 def __init__(self, geometry: shapely.Polygon | shapely.MultiPolygon) -> None:
44 self._impl = geometry
46 @property
47 def area(self) -> float:
48 """The area of the region (`float`)."""
49 return self._impl.area
51 @property
52 def bbox(self) -> Box:
53 """The integer-coordinate bounding box of the region (`Box`).
55 Because a `Box` logically contains the entirety of the pixels on its
56 boundary, but the centers of those pixels are the numerical values of
57 its bounds, the region may contain vertices that are up to 0.5 beyond
58 the integer box coordinates in either dimension.
59 """
60 x_min, y_min, x_max, y_max = self._impl.bounds
61 return Box.factory[
62 round_half_up(y_min) : round_half_down(y_max) + 1,
63 round_half_up(x_min) : round_half_down(x_max) + 1,
64 ]
66 @property
67 def wkt(self) -> str:
68 """The 'Well-Known Text' representation of this region (`str`)."""
69 return self._impl.wkt
71 @staticmethod
72 def from_wkt(wkt: str) -> Region:
73 """Construct from a 'Well-Known Text' string.
75 Parameters
76 ----------
77 wkt
78 Well-Known Text representation of the region.
79 """
80 impl = shapely.from_wkt(wkt)
81 if not isinstance(impl, shapely.Polygon | shapely.MultiPolygon): 81 ↛ 82line 81 didn't jump to line 82 because the condition on line 81 was never true
82 raise ValueError("Only Polygon and MultiPolygon geometries can be converted to Regions.")
83 return Region(impl).try_to_polygon()
85 def __str__(self) -> str:
86 return self._impl.wkt
88 def __repr__(self) -> str:
89 return f"Region.from_wkt({self._impl.wkt!r})"
91 def __eq__(self, other: object) -> bool:
92 if isinstance(other, Region):
93 return bool(shapely.equals(self._impl, other._impl))
94 return NotImplemented
96 @overload
97 def contains(self, other: Polygon) -> bool: ... 97 ↛ exitline 97 didn't return from function 'contains' because
99 @overload
100 def contains(self, *, x: int, y: int) -> bool: ... 100 ↛ exitline 100 didn't return from function 'contains' because
102 @overload
103 def contains(self, *, x: float, y: float) -> bool: ... 103 ↛ exitline 103 didn't return from function 'contains' because
105 @overload
106 def contains(self, *, x: np.ndarray, y: np.ndarray) -> np.ndarray: ... 106 ↛ exitline 106 didn't return from function 'contains' because
108 def contains(
109 self,
110 other: Region | None = None,
111 *,
112 x: float | int | np.ndarray | None = None,
113 y: float | int | np.ndarray | None = None,
114 ) -> bool | np.ndarray:
115 """Test whether the geometry contains the given points or another
116 geometry.
118 Parameters
119 ----------
120 other
121 Another geometry to compare to. Not compatible with the ``y`` and
122 ``x`` arguments.
123 x
124 One or more floating-point or integer X coordinates to test for
125 containment. If an array, an array of results will be returned.
126 y
127 One or more floating-point or integer Y coordinates to test for
128 containment. If an array, an array of results will be returned.
129 """
130 if other is not None:
131 if x is not None or y is not None: 131 ↛ 132line 131 didn't jump to line 132 because the condition on line 131 was never true
132 raise TypeError("Too many arguments to 'Region.contains'.")
133 return self._impl.contains(other._impl)
134 elif x is None or y is None: 134 ↛ 135line 134 didn't jump to line 135 because the condition on line 134 was never true
135 raise TypeError("Not enough arguments to 'Region.contains'.")
136 else:
137 # Quibbles about bool vs numpy.bool_ as the return type.
138 return shapely.contains_xy(self._impl, x=x, y=y) # type: ignore[return-value]
140 def intersection(self, other: Bounds) -> Bounds:
141 """Compute the intersection of this region with a `Bounds` object.
143 Notes
144 -----
145 Because `Region` implements the `Bounds` interface, its intersections
146 need to support all other `Bounds` objects. This is not true of other
147 `Region` point-set operations like `union` and `difference`.
149 Parameters
150 ----------
151 other
152 Bounds to intersect with this region.
153 """
154 from ._concrete_bounds import _intersect_region
156 return _intersect_region(self, other)
158 def union(self, other: Region) -> Region:
159 """Compute the point-set union of this region with another.
161 Parameters
162 ----------
163 other
164 Region to union with this one.
165 """
166 impl = shapely.union(self._impl, other._impl)
167 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
168 "A union of Polygons and MultiPolygons should be one of those."
169 )
170 return Region(impl).try_to_polygon()
172 def difference(self, other: Region) -> Region:
173 """Compute the point-set difference of this region with another.
175 Parameters
176 ----------
177 other
178 Region to subtract from this one.
179 """
180 impl = shapely.difference(self._impl, other._impl)
181 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
182 "A difference of Polygons and MultiPolygons should be one of those."
183 )
184 return Region(impl).try_to_polygon()
186 def try_to_polygon(self) -> Region:
187 """If the underlying geometry is a single polygon with no holes,
188 return a `Polygon` instance holding it.
190 In all other cases ``self`` is returned.
191 """
192 impl = self._impl
193 if isinstance(impl, shapely.MultiPolygon) and len(impl.geoms) == 1: 193 ↛ 194line 193 didn't jump to line 194 because the condition on line 193 was never true
194 impl = impl.geoms[0]
195 if isinstance(impl, shapely.Polygon) and not impl.interiors:
196 vertices = np.array(impl.exterior.coords)
197 return Polygon(x_vertices=vertices[:, 0], y_vertices=vertices[:, 1])
198 return self
200 def try_to_box(self) -> Region | Box:
201 """If the underlying geometry is a rectangle that fully covers integer
202 pixels (i.e. has all vertices at half-integer positions), return the
203 equivalent `Box`.
205 In all other cases ``self`` is returned.
206 """
207 if Polygon.from_box(self.bbox) == self:
208 return self.bbox
209 return self
211 def to_shapely(self) -> shapely.Polygon | shapely.MultiPolygon:
212 """Convert to a `shapely.Polygon` or `shapely.MultiPolygon` object."""
213 return self._impl
215 def serialize(self) -> RegionSerializationModel:
216 """Serialize the region to a Pydantic model.
218 Region serialization uses a subset of the GeoJSON specification (IETF
219 RFC 7946).
220 """
221 return RegionSerializationModel.model_validate_json(shapely.to_geojson(self._impl))
224class Polygon(Region):
225 """A simple 2-d polygon in Euclidean coordinates, with no holes.
227 Parameters
228 ----------
229 x_vertices
230 The x coordinates of the vertices of the polygon.
231 y_vertices
232 The y coordinate of the vertices of the polygon.
233 """
235 def __init__(self, *, x_vertices: npt.ArrayLike, y_vertices: npt.ArrayLike) -> None:
236 self._vertices = np.stack(
237 [np.asarray(x_vertices).flat, np.asarray(y_vertices).flat], dtype=np.float64
238 ).transpose()
239 self._vertices.flags.writeable = False
240 super().__init__(shapely.Polygon(self._vertices))
242 @staticmethod
243 def from_box(box: Box) -> Polygon:
244 """Construct from an integer-coordinate box.
246 Parameters
247 ----------
248 box
249 Integer-coordinate box to convert to a polygon.
251 Notes
252 -----
253 Because the integer min and max coordinates of the box are
254 interpreted as pixel centers, these are expanded by 0.5 on all sides
255 before using them to form the polygon vertices.
256 """
257 return Polygon(
258 x_vertices=[box.x.min - 0.5, box.x.min - 0.5, box.x.max + 0.5, box.x.max + 0.5],
259 y_vertices=[box.y.min - 0.5, box.y.max + 0.5, box.y.max + 0.5, box.y.min - 0.5],
260 )
262 @property
263 def n_vertices(self) -> int:
264 """The number of vertices in the polygon."""
265 return self._vertices.shape[0]
267 @property
268 def x_vertices(self) -> np.ndarray:
269 """The x coordinates of the vertices of the polygon.
271 This is a read-only array; polygons are immutable.
272 """
273 return self._vertices[:, 0]
275 @property
276 def y_vertices(self) -> np.ndarray:
277 """The y coordinates of the vertices of the polygon.
279 This is a read-only array; polygons are immutable.
280 """
281 return self._vertices[:, 1]
283 def __repr__(self) -> str:
284 return f"Polygon(x_vertices={self.x_vertices!r}, y_vertices={self.y_vertices!r})"
286 @staticmethod
287 def from_legacy(legacy: LegacyPolygon) -> Polygon:
288 """Convert from a legacy `lsst.afw.geom.Polygon` instance.
290 Parameters
291 ----------
292 legacy
293 Legacy `lsst.afw.geom.Polygon` to convert.
294 """
295 vertices = legacy.getVertices()
296 x_vertices = np.zeros(len(vertices), dtype=np.float64)
297 y_vertices = np.zeros(len(vertices), dtype=np.float64)
298 for n, point in enumerate(vertices):
299 x_vertices[n] = point.x
300 y_vertices[n] = point.y
301 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices)
303 def to_legacy(self) -> LegacyPolygon:
304 """Convert to a legacy `lsst.afw.geom.Polygon` instance."""
305 from lsst.afw.geom import Polygon as LegacyPolygon
306 from lsst.geom import Point2D
308 return LegacyPolygon([Point2D(x, y) for x, y in zip(self.x_vertices, self.y_vertices)])
311class RegionSerializationModel(pydantic.BaseModel):
312 """Serialization model for `Region` and `Polygon`.
314 This model is a subset of the GeoJSON specification (IETF RFC 7946).
315 """
317 type: Literal["Polygon", "MultiPolygon"] = pydantic.Field(description="Geometry type.")
319 coordinates: list[list[tuple[float, float] | list[tuple[float, float]]]] = pydantic.Field(
320 description="Vertices of the polygon or polygons."
321 )
323 def deserialize(self) -> Region:
324 """Deserialize into a `Region` (a `Polygon`, if possible)."""
325 region_impl = shapely.from_geojson(self.model_dump_json())
326 assert isinstance(region_impl, shapely.Polygon | shapely.MultiPolygon), (
327 "Other geometry types are not used."
328 )
329 return Region(region_impl).try_to_polygon()