Coverage for python/lsst/images/_polygon.py: 42%
118 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-05-30 09:08 +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):
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."""
74 impl = shapely.from_wkt(wkt)
75 if not isinstance(impl, shapely.Polygon | shapely.MultiPolygon):
76 raise ValueError("Only Polygon and MultiPolygon geometries can be converted to Regions.")
77 return Region(impl).try_to_polygon()
79 def __str__(self) -> str:
80 return self._impl.wkt
82 def __repr__(self) -> str:
83 return f"Region.from_wkt({self._impl.wkt!r})"
85 def __eq__(self, other: object) -> bool:
86 if isinstance(other, Region):
87 return bool(shapely.equals(self._impl, other._impl))
88 return NotImplemented
90 @overload
91 def contains(self, other: Polygon) -> bool: ... 91 ↛ exitline 91 didn't return from function 'contains' because
93 @overload
94 def contains(self, *, x: int, y: int) -> bool: ... 94 ↛ exitline 94 didn't return from function 'contains' because
96 @overload
97 def contains(self, *, x: float, y: float) -> bool: ... 97 ↛ exitline 97 didn't return from function 'contains' because
99 @overload
100 def contains(self, *, x: np.ndarray, y: np.ndarray) -> np.ndarray: ... 100 ↛ exitline 100 didn't return from function 'contains' because
102 def contains(
103 self,
104 other: Region | None = None,
105 *,
106 x: float | int | np.ndarray | None = None,
107 y: float | int | np.ndarray | None = None,
108 ) -> bool | np.ndarray:
109 """Test whether the geometry contains the given points or another
110 geometry.
112 Parameters
113 ----------
114 other
115 Another geometry to compare to. Not compatible with the ``y`` and
116 ``x`` arguments.
117 x
118 One or more floating-point or integer X coordinates to test for
119 containment. If an array, an array of results will be returned.
120 y
121 One or more floating-point or integer Y coordinates to test for
122 containment. If an array, an array of results will be returned.
123 """
124 if other is not None:
125 if x is not None or y is not None:
126 raise TypeError("Too many arguments to 'Region.contains'.")
127 return self._impl.contains(other._impl)
128 elif x is None or y is None:
129 raise TypeError("Not enough arguments to 'Region.contains'.")
130 else:
131 # Quibbles about bool vs numpy.bool_ as the return type.
132 return shapely.contains_xy(self._impl, x=x, y=y) # type: ignore[return-value]
134 def intersection(self, other: Bounds) -> Bounds:
135 """Compute the intersection of this region with a `Bounds` object.
137 Notes
138 -----
139 Because `Region` implements the `Bounds` interface, its intersections
140 need to support all other `Bounds` objects. This is not true of other
141 `Region` point-set operations like `union` and `difference`.
142 """
143 from ._concrete_bounds import _intersect_region
145 return _intersect_region(self, other)
147 def union(self, other: Region) -> Region:
148 """Compute the point-set union of this region with another."""
149 impl = shapely.union(self._impl, other._impl)
150 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
151 "A union of Polygons and MultiPolygons should be one of those."
152 )
153 return Region(impl).try_to_polygon()
155 def difference(self, other: Region) -> Region:
156 """Compute the point-set difference of this region with another."""
157 impl = shapely.difference(self._impl, other._impl)
158 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
159 "A difference of Polygons and MultiPolygons should be one of those."
160 )
161 return Region(impl).try_to_polygon()
163 def try_to_polygon(self) -> Region:
164 """If the underlying geometry is a single polygon with no holes,
165 return a `Polygon` instance holding it.
167 In all other cases ``self`` is returned.
168 """
169 impl = self._impl
170 if isinstance(impl, shapely.MultiPolygon) and len(impl.geoms) == 1:
171 impl = impl.geoms[0]
172 if isinstance(impl, shapely.Polygon) and not impl.interiors:
173 vertices = np.array(impl.exterior.coords)
174 return Polygon(x_vertices=vertices[:, 0], y_vertices=vertices[:, 1])
175 return self
177 def try_to_box(self) -> Region | Box:
178 """If the underlying geometry is a rectangle that fully covers integer
179 pixels (i.e. has all vertices at half-integer positions), return the
180 equivalent `Box`.
182 In all other cases ``self`` is returned.
183 """
184 if Polygon.from_box(self.bbox) == self:
185 return self.bbox
186 return self
188 def to_shapely(self) -> shapely.Polygon | shapely.MultiPolygon:
189 """Convert to a `shapely.Polygon` or `shapely.MultiPolygon` object."""
190 return self._impl
192 def serialize(self) -> RegionSerializationModel:
193 """Serialize the region to a Pydantic model.
195 Region serialization uses a subset of the GeoJSON specification (IETF
196 RFC 7946).
197 """
198 return RegionSerializationModel.model_validate_json(shapely.to_geojson(self._impl))
201class Polygon(Region):
202 """A simple 2-d polygon in Euclidean coordinates, with no holes.
204 Parameters
205 ----------
206 x_vertices
207 The x coordinates of the vertices of the polygon.
208 y_vertices
209 The y coordinate of the vertices of the polygon.
210 """
212 def __init__(self, *, x_vertices: npt.ArrayLike, y_vertices: npt.ArrayLike):
213 self._vertices = np.stack(
214 [np.asarray(x_vertices).flat, np.asarray(y_vertices).flat], dtype=np.float64
215 ).transpose()
216 self._vertices.flags.writeable = False
217 super().__init__(shapely.Polygon(self._vertices))
219 @staticmethod
220 def from_box(box: Box) -> Polygon:
221 """Construct from an integer-coordinate box.
223 Notes
224 -----
225 Because the integer min and max coordinates of the box are
226 interpreted as pixel centers, these are expanded by 0.5 on all sides
227 before using them to form the polygon vertices.
228 """
229 return Polygon(
230 x_vertices=[box.x.min - 0.5, box.x.min - 0.5, box.x.max + 0.5, box.x.max + 0.5],
231 y_vertices=[box.y.min - 0.5, box.y.max + 0.5, box.y.max + 0.5, box.y.min - 0.5],
232 )
234 @property
235 def n_vertices(self) -> int:
236 """The number of vertices in the polygon."""
237 return self._vertices.shape[0]
239 @property
240 def x_vertices(self) -> np.ndarray:
241 """The x coordinates of the vertices of the polygon.
243 This is a read-only array; polygons are immutable.
244 """
245 return self._vertices[:, 0]
247 @property
248 def y_vertices(self) -> np.ndarray:
249 """The y coordinates of the vertices of the polygon.
251 This is a read-only array; polygons are immutable.
252 """
253 return self._vertices[:, 1]
255 def __repr__(self) -> str:
256 return f"Polygon(x_vertices={self.x_vertices!r}, y_vertices={self.y_vertices!r})"
258 @staticmethod
259 def from_legacy(legacy: LegacyPolygon) -> Polygon:
260 """Convert from a legacy `lsst.afw.geom.Polygon` instance."""
261 vertices = legacy.getVertices()
262 x_vertices = np.zeros(len(vertices), dtype=np.float64)
263 y_vertices = np.zeros(len(vertices), dtype=np.float64)
264 for n, point in enumerate(vertices):
265 x_vertices[n] = point.x
266 y_vertices[n] = point.y
267 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices)
269 def to_legacy(self) -> LegacyPolygon:
270 """Convert to a legacy `lsst.afw.geom.Polygon` instance."""
271 from lsst.afw.geom import Polygon as LegacyPolygon
272 from lsst.geom import Point2D
274 return LegacyPolygon([Point2D(x, y) for x, y in zip(self.x_vertices, self.y_vertices)])
277class RegionSerializationModel(pydantic.BaseModel):
278 """Serialization model for `Region` and `Polygon`.
280 This model is a subset of the GeoJSON specification (IETF RFC 7946).
281 """
283 type: Literal["Polygon", "MultiPolygon"] = pydantic.Field(description="Geometry type.")
285 coordinates: list[list[tuple[float, float] | list[tuple[float, float]]]] = pydantic.Field(
286 description="Vertices of the polygon or polygons."
287 )
289 def deserialize(self) -> Region:
290 """Deserialize into a `Region` (a `Polygon`, if possible)."""
291 region_impl = shapely.from_geojson(self.model_dump_json())
292 assert isinstance(region_impl, shapely.Polygon | shapely.MultiPolygon), (
293 "Other geometry types are not used."
294 )
295 return Region(region_impl).try_to_polygon()