Coverage for python/lsst/images/_polygon.py: 91%
145 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-07-09 02:27 -0700
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 pydantic_core.core_schema as pcs
22import shapely
23from pydantic.json_schema import JsonSchemaValue
25from ._geom import XY, Bounds, Box
26from .utils import round_half_down, round_half_up
28if TYPE_CHECKING:
29 from ._transforms import Transform
31 try:
32 from lsst.afw.geom import Polygon as LegacyPolygon
33 except ImportError:
34 type LegacyPolygon = Any # type: ignore[no-redef]
37class Region:
38 """A 2-d Euclidean region represented as one or more polygons with
39 optional holes.
41 Parameters
42 ----------
43 geometry
44 A polygon or multi-polygon from the Shapely library.
45 """
47 def __init__(self, geometry: shapely.Polygon | shapely.MultiPolygon) -> None:
48 self._impl = geometry
50 @property
51 def area(self) -> float:
52 """The area of the region (`float`)."""
53 return self._impl.area
55 @property
56 def bbox(self) -> Box:
57 """The integer-coordinate bounding box of the region (`Box`).
59 Because a `Box` logically contains the entirety of the pixels on its
60 boundary, but the centers of those pixels are the numerical values of
61 its bounds, the region may contain vertices that are up to 0.5 beyond
62 the integer box coordinates in either dimension.
63 """
64 x_min, y_min, x_max, y_max = self._impl.bounds
65 return Box.factory[
66 round_half_up(y_min) : round_half_down(y_max) + 1,
67 round_half_up(x_min) : round_half_down(x_max) + 1,
68 ]
70 @property
71 def wkt(self) -> str:
72 """The 'Well-Known Text' representation of this region (`str`)."""
73 return self._impl.wkt
75 @staticmethod
76 def from_wkt(wkt: str) -> Region:
77 """Construct from a 'Well-Known Text' string.
79 Parameters
80 ----------
81 wkt
82 Well-Known Text representation of the region.
83 """
84 impl = shapely.from_wkt(wkt)
85 if not isinstance(impl, shapely.Polygon | shapely.MultiPolygon): 85 ↛ 86line 85 didn't jump to line 86 because the condition on line 85 was never true
86 raise ValueError("Only Polygon and MultiPolygon geometries can be converted to Regions.")
87 return Region(impl).try_to_polygon()
89 def __str__(self) -> str:
90 return self._impl.wkt
92 def __repr__(self) -> str:
93 return f"Region.from_wkt({self._impl.wkt!r})"
95 def __eq__(self, other: object) -> bool:
96 if isinstance(other, Region):
97 return bool(shapely.equals(self._impl, other._impl))
98 return NotImplemented
100 @overload
101 def contains(self, other: Polygon) -> bool: ... 101 ↛ exitline 101 didn't return from function 'contains' because
103 @overload
104 def contains(self, *, x: int, y: int) -> bool: ... 104 ↛ exitline 104 didn't return from function 'contains' because
106 @overload
107 def contains(self, *, x: float, y: float) -> bool: ... 107 ↛ exitline 107 didn't return from function 'contains' because
109 @overload
110 def contains(self, *, x: np.ndarray, y: np.ndarray) -> np.ndarray: ... 110 ↛ exitline 110 didn't return from function 'contains' because
112 def contains(
113 self,
114 other: Region | None = None,
115 *,
116 x: float | int | np.ndarray | None = None,
117 y: float | int | np.ndarray | None = None,
118 ) -> bool | np.ndarray:
119 """Test whether the geometry contains the given points or another
120 geometry.
122 Parameters
123 ----------
124 other
125 Another geometry to compare to. Not compatible with the ``y`` and
126 ``x`` arguments.
127 x
128 One or more floating-point or integer X coordinates to test for
129 containment. If an array, an array of results will be returned.
130 y
131 One or more floating-point or integer Y coordinates to test for
132 containment. If an array, an array of results will be returned.
133 """
134 if other is not None:
135 if x is not None or y is not None: 135 ↛ 136line 135 didn't jump to line 136 because the condition on line 135 was never true
136 raise TypeError("Too many arguments to 'Region.contains'.")
137 return self._impl.contains(other._impl)
138 elif x is None or y is None: 138 ↛ 139line 138 didn't jump to line 139 because the condition on line 138 was never true
139 raise TypeError("Not enough arguments to 'Region.contains'.")
140 else:
141 # Quibbles about bool vs numpy.bool_ as the return type.
142 return shapely.contains_xy(self._impl, x=x, y=y) # type: ignore[return-value]
144 @overload
145 def intersection(self, other: Region) -> Region: ... 145 ↛ exitline 145 didn't return from function 'intersection' because
147 @overload
148 def intersection(self, other: Box) -> Region | Box: ... 148 ↛ exitline 148 didn't return from function 'intersection' because
150 @overload
151 def intersection(self, other: Bounds) -> Bounds: ... 151 ↛ exitline 151 didn't return from function 'intersection' because
153 def intersection(self, other: Bounds) -> Bounds:
154 """Compute the intersection of this region with a `Bounds` object.
156 Notes
157 -----
158 Because `Region` implements the `Bounds` interface, its intersections
159 need to support all other `Bounds` objects. This is not true of other
160 `Region` point-set operations like `union` and `difference`.
162 Parameters
163 ----------
164 other
165 Bounds to intersect with this region.
166 """
167 from ._concrete_bounds import _intersect_region
169 return _intersect_region(self, other)
171 def union(self, other: Region) -> Region:
172 """Compute the point-set union of this region with another.
174 Parameters
175 ----------
176 other
177 Region to union with this one.
178 """
179 impl = shapely.union(self._impl, other._impl)
180 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
181 "A union of Polygons and MultiPolygons should be one of those."
182 )
183 return Region(impl).try_to_polygon()
185 def difference(self, other: Region) -> Region:
186 """Compute the point-set difference of this region with another.
188 Parameters
189 ----------
190 other
191 Region to subtract from this one.
192 """
193 impl = shapely.difference(self._impl, other._impl)
194 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
195 "A difference of Polygons and MultiPolygons should be one of those."
196 )
197 return Region(impl).try_to_polygon()
199 def try_to_polygon(self) -> Region:
200 """If the underlying geometry is a single polygon with no holes,
201 return a `Polygon` instance holding it.
203 In all other cases ``self`` is returned.
204 """
205 impl = self._impl
206 if isinstance(impl, shapely.MultiPolygon) and len(impl.geoms) == 1: 206 ↛ 207line 206 didn't jump to line 207 because the condition on line 206 was never true
207 impl = impl.geoms[0]
208 if isinstance(impl, shapely.Polygon) and not impl.interiors:
209 vertices = np.array(impl.exterior.coords)
210 return Polygon(x_vertices=vertices[:, 0], y_vertices=vertices[:, 1])
211 return self
213 def try_to_box(self) -> Region | Box:
214 """If the underlying geometry is a rectangle that fully covers integer
215 pixels (i.e. has all vertices at half-integer positions), return the
216 equivalent `Box`.
218 In all other cases ``self`` is returned.
219 """
220 if Polygon.from_box(self.bbox) == self:
221 return self.bbox
222 return self
224 def to_shapely(self) -> shapely.Polygon | shapely.MultiPolygon:
225 """Convert to a `shapely.Polygon` or `shapely.MultiPolygon` object."""
226 return self._impl
228 def transform(self, transform: Transform[Any, Any]) -> Region:
229 """Transform the coordinates of the region, returning a new one.
231 Parameters
232 ----------
233 transform
234 Coordinate transform to apply (in the forward direction).
236 Notes
237 -----
238 This applies the transform to all vertices, assuming that the
239 transform is close enough to affine that the topology of the geometry
240 does not change and straight-line edges do not need to be subsampled.
241 """
243 def wrapper(x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:
244 return transform.apply_forward(x=x, y=y)
246 return Region(
247 # Shapely overloads don't seem to have been annotated rigorously
248 shapely.transform(self._impl, wrapper, interleaved=False) # type: ignore[arg-type]
249 ).try_to_polygon()
251 def serialize(self) -> RegionSerializationModel:
252 """Serialize the region to a Pydantic model.
254 Region serialization uses a subset of the GeoJSON specification (IETF
255 RFC 7946).
256 """
257 return RegionSerializationModel.model_validate_json(shapely.to_geojson(self._impl))
259 @classmethod
260 def __get_pydantic_core_schema__(
261 cls, source_type: Any, handler: pydantic.GetCoreSchemaHandler
262 ) -> pcs.CoreSchema:
263 from_model_schema = pcs.chain_schema(
264 [
265 handler(RegionSerializationModel),
266 pcs.no_info_plain_validator_function(RegionSerializationModel.deserialize),
267 ]
268 )
269 return pcs.json_or_python_schema(
270 json_schema=from_model_schema,
271 python_schema=pcs.union_schema([pcs.is_instance_schema(cls), from_model_schema]),
272 serialization=pcs.plain_serializer_function_ser_schema(cls.serialize),
273 )
275 @classmethod
276 def __get_pydantic_json_schema__(
277 cls, schema: pcs.CoreSchema, handler: pydantic.GetJsonSchemaHandler
278 ) -> JsonSchemaValue:
279 return handler(RegionSerializationModel.__pydantic_core_schema__)
282class Polygon(Region):
283 """A simple 2-d polygon in Euclidean coordinates, with no holes.
285 Parameters
286 ----------
287 x_vertices
288 The x coordinates of the vertices of the polygon.
289 y_vertices
290 The y coordinate of the vertices of the polygon.
291 """
293 def __init__(self, *, x_vertices: npt.ArrayLike, y_vertices: npt.ArrayLike) -> None:
294 self._vertices = np.stack(
295 [np.asarray(x_vertices).flat, np.asarray(y_vertices).flat], dtype=np.float64
296 ).transpose()
297 self._vertices.flags.writeable = False
298 super().__init__(shapely.Polygon(self._vertices))
300 @staticmethod
301 def from_box(box: Box) -> Polygon:
302 """Construct from an integer-coordinate box.
304 Parameters
305 ----------
306 box
307 Integer-coordinate box to convert to a polygon.
309 Notes
310 -----
311 Because the integer min and max coordinates of the box are
312 interpreted as pixel centers, these are expanded by 0.5 on all sides
313 before using them to form the polygon vertices.
314 """
315 return Polygon(
316 x_vertices=[box.x.min - 0.5, box.x.min - 0.5, box.x.max + 0.5, box.x.max + 0.5],
317 y_vertices=[box.y.min - 0.5, box.y.max + 0.5, box.y.max + 0.5, box.y.min - 0.5],
318 )
320 @property
321 def n_vertices(self) -> int:
322 """The number of vertices in the polygon."""
323 return self._vertices.shape[0]
325 @property
326 def x_vertices(self) -> np.ndarray:
327 """The x coordinates of the vertices of the polygon.
329 This is a read-only array; polygons are immutable.
330 """
331 return self._vertices[:, 0]
333 @property
334 def y_vertices(self) -> np.ndarray:
335 """The y coordinates of the vertices of the polygon.
337 This is a read-only array; polygons are immutable.
338 """
339 return self._vertices[:, 1]
341 @property
342 def centroid(self) -> XY[float]:
343 """The centroid of the polygon (`XY` [`float`])."""
344 c = self._impl.centroid
345 return XY(x=c.x, y=c.y)
347 def __repr__(self) -> str:
348 return f"Polygon(x_vertices={self.x_vertices!r}, y_vertices={self.y_vertices!r})"
350 def transform(self, transform: Transform[Any, Any]) -> Polygon:
351 # Docstring inherited.
352 result = super().transform(transform)
353 assert isinstance(result, Polygon), "Transforming a polygon should not change its topology."
354 return result
356 @staticmethod
357 def from_legacy(legacy: LegacyPolygon) -> Polygon:
358 """Convert from a legacy `lsst.afw.geom.Polygon` instance.
360 Parameters
361 ----------
362 legacy
363 Legacy `lsst.afw.geom.Polygon` to convert.
364 """
365 vertices = legacy.getVertices()
366 x_vertices = np.zeros(len(vertices), dtype=np.float64)
367 y_vertices = np.zeros(len(vertices), dtype=np.float64)
368 for n, point in enumerate(vertices):
369 x_vertices[n] = point.x
370 y_vertices[n] = point.y
371 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices)
373 def to_legacy(self) -> LegacyPolygon:
374 """Convert to a legacy `lsst.afw.geom.Polygon` instance."""
375 from lsst.afw.geom import Polygon as LegacyPolygon
376 from lsst.geom import Point2D
378 return LegacyPolygon([Point2D(x, y) for x, y in zip(self.x_vertices, self.y_vertices)])
381class RegionSerializationModel(pydantic.BaseModel):
382 """Serialization model for `Region` and `Polygon`.
384 This model is a subset of the GeoJSON specification (IETF RFC 7946).
385 """
387 type: Literal["Polygon", "MultiPolygon"] = pydantic.Field(description="Geometry type.")
389 coordinates: list[list[tuple[float, float] | list[tuple[float, float]]]] = pydantic.Field(
390 description="Vertices of the polygon or polygons."
391 )
393 def deserialize(self) -> Region:
394 """Deserialize into a `Region` (a `Polygon`, if possible)."""
395 region_impl = shapely.from_geojson(self.model_dump_json())
396 assert isinstance(region_impl, shapely.Polygon | shapely.MultiPolygon), (
397 "Other geometry types are not used."
398 )
399 return Region(region_impl).try_to_polygon()