Coverage for python/lsst/images/_polygon.py: 47%
158 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-30 08:38 +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, assert_type, 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, YX, Bounds, Box
26from .describe import DescribableMixin, DescribeOptions, Report, ReportField
27from .utils import round_half_down, round_half_up
29if TYPE_CHECKING:
30 from ._geom import Bounds
31 from ._transforms import Transform
33 try:
34 from lsst.afw.geom import Polygon as LegacyPolygon
35 except ImportError:
36 type LegacyPolygon = Any # type: ignore[no-redef]
39class Region(DescribableMixin):
40 """A 2-d Euclidean region represented as one or more polygons with
41 optional holes.
43 Parameters
44 ----------
45 geometry
46 A polygon or multi-polygon from the Shapely library.
47 """
49 def __init__(self, geometry: shapely.Polygon | shapely.MultiPolygon) -> None:
50 self._impl = geometry
52 @property
53 def area(self) -> float:
54 """The area of the region (`float`)."""
55 return self._impl.area
57 @property
58 def bbox(self) -> Box:
59 """The integer-coordinate bounding box of the region (`Box`).
61 Because a `Box` logically contains the entirety of the pixels on its
62 boundary, but the centers of those pixels are the numerical values of
63 its bounds, the region may contain vertices that are up to 0.5 beyond
64 the integer box coordinates in either dimension.
65 """
66 x_min, y_min, x_max, y_max = self._impl.bounds
67 return Box.factory[
68 round_half_up(y_min) : round_half_down(y_max) + 1,
69 round_half_up(x_min) : round_half_down(x_max) + 1,
70 ]
72 @property
73 def wkt(self) -> str:
74 """The 'Well-Known Text' representation of this region (`str`)."""
75 return self._impl.wkt
77 @staticmethod
78 def from_wkt(wkt: str) -> Region:
79 """Construct from a 'Well-Known Text' string.
81 Parameters
82 ----------
83 wkt
84 Well-Known Text representation of the region.
85 """
86 impl = shapely.from_wkt(wkt)
87 if not isinstance(impl, shapely.Polygon | shapely.MultiPolygon): 87 ↛ 88line 87 didn't jump to line 88 because the condition on line 87 was never true
88 raise ValueError("Only Polygon and MultiPolygon geometries can be converted to Regions.")
89 return Region(impl).try_to_polygon()
91 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:
92 """Return a `Report` describing this region.
94 Parameters
95 ----------
96 options : `DescribeOptions`, optional
97 Unused; accepted for interface compatibility.
98 """
99 return Report(
100 type_name="Region",
101 summary=self._impl.wkt,
102 fields=[
103 ReportField(
104 label="wkt",
105 value=self._impl.wkt,
106 repr_value=repr(self._impl.wkt),
107 positional=True,
108 )
109 ],
110 )
112 def __repr__(self) -> str:
113 return f"Region.from_wkt({self._impl.wkt!r})"
115 def __eq__(self, other: object) -> bool:
116 if isinstance(other, Region):
117 return bool(shapely.equals(self._impl, other._impl))
118 return NotImplemented
120 @overload
121 def contains(self, point: XY[int | float] | YX[int | float], /) -> bool: ...
123 @overload
124 def contains(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> np.ndarray: ...
126 @overload
127 def contains(self, /, *, x: int | float, y: int | float) -> bool: ...
129 @overload
130 def contains(self, /, *, x: npt.ArrayLike, y: npt.ArrayLike) -> np.ndarray: ...
132 @overload
133 def contains(self, other: Polygon, /) -> bool: ...
135 def contains(
136 self,
137 other: Region | XY[Any] | YX[Any] | None = None,
138 /,
139 *,
140 x: Any = None,
141 y: Any = None,
142 ) -> bool | np.ndarray:
143 """Test whether the geometry contains the given points or another
144 geometry.
146 Parameters
147 ----------
148 other
149 Another geometry, or an `XY` or `YX` coordinate pair, to compare
150 to. Mutually exclusive with ``x`` and ``y``.
151 x
152 One or more X coordinates to test for containment, as a scalar or
153 any array-like. Results are broadcast against ``y``.
154 Mutually exclusive with ``other``.
155 y
156 One or more Y coordinates to test for containment, as a scalar or
157 any array-like. Results are broadcast against ``x``.
158 Mutually exclusive with ``other``.
159 """
160 match other:
161 case None:
162 if x is None or y is None: 162 ↛ 163line 162 didn't jump to line 163 because the condition on line 162 was never true
163 raise TypeError("Pass either a point or both x= and y= to 'Region.contains'.")
164 case Region():
165 if x is not None or y is not None: 165 ↛ 166line 165 didn't jump to line 166 because the condition on line 165 was never true
166 raise TypeError("'Region.contains' point argument is mutually exclusive with x= and y=.")
167 return self._impl.contains(other._impl)
168 case XY() | YX(): 168 ↛ 172line 168 didn't jump to line 172 because the pattern on line 168 always matched
169 if x is not None or y is not None:
170 raise TypeError("'Region.contains' point argument is mutually exclusive with x= and y=.")
171 x, y = other.x, other.y
172 case _:
173 raise TypeError(f"Unexpected positional argument type: {type(other)!r}.")
174 return shapely.contains_xy(self._impl, x=x, y=y)
176 @overload
177 def intersection(self, other: Region) -> Region: ...
179 @overload
180 def intersection(self, other: Box) -> Region | Box: ...
182 @overload
183 def intersection(self, other: Bounds) -> Bounds: ...
185 def intersection(self, other: Bounds) -> Bounds:
186 """Compute the intersection of this region with a `Bounds` object.
188 Notes
189 -----
190 Because `Region` implements the `Bounds` interface, its intersections
191 need to support all other `Bounds` objects. This is not true of other
192 `Region` point-set operations like `union` and `difference`.
194 Parameters
195 ----------
196 other
197 Bounds to intersect with this region.
198 """
199 from ._concrete_bounds import _intersect_region
201 return _intersect_region(self, other)
203 def union(self, other: Region) -> Region:
204 """Compute the point-set union of this region with another.
206 Parameters
207 ----------
208 other
209 Region to union with this one.
210 """
211 impl = shapely.union(self._impl, other._impl)
212 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
213 "A union of Polygons and MultiPolygons should be one of those."
214 )
215 return Region(impl).try_to_polygon()
217 def difference(self, other: Region) -> Region:
218 """Compute the point-set difference of this region with another.
220 Parameters
221 ----------
222 other
223 Region to subtract from this one.
224 """
225 impl = shapely.difference(self._impl, other._impl)
226 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), (
227 "A difference of Polygons and MultiPolygons should be one of those."
228 )
229 return Region(impl).try_to_polygon()
231 def try_to_polygon(self) -> Region:
232 """If the underlying geometry is a single polygon with no holes,
233 return a `Polygon` instance holding it.
235 In all other cases ``self`` is returned.
236 """
237 impl = self._impl
238 if isinstance(impl, shapely.MultiPolygon) and len(impl.geoms) == 1: 238 ↛ 239line 238 didn't jump to line 239 because the condition on line 238 was never true
239 impl = impl.geoms[0]
240 if isinstance(impl, shapely.Polygon) and not impl.interiors:
241 vertices = np.array(impl.exterior.coords)
242 return Polygon(x_vertices=vertices[:, 0], y_vertices=vertices[:, 1])
243 return self
245 def try_to_box(self) -> Region | Box:
246 """If the underlying geometry is a rectangle that fully covers integer
247 pixels (i.e. has all vertices at half-integer positions), return the
248 equivalent `Box`.
250 In all other cases ``self`` is returned.
251 """
252 if Polygon.from_box(self.bbox) == self:
253 return self.bbox
254 return self
256 def to_shapely(self) -> shapely.Polygon | shapely.MultiPolygon:
257 """Convert to a `shapely.Polygon` or `shapely.MultiPolygon` object."""
258 return self._impl
260 def transform(self, transform: Transform[Any, Any]) -> Region:
261 """Transform the coordinates of the region, returning a new one.
263 Parameters
264 ----------
265 transform
266 Coordinate transform to apply (in the forward direction).
268 Notes
269 -----
270 This applies the transform to all vertices, assuming that the
271 transform is close enough to affine that the topology of the geometry
272 does not change and straight-line edges do not need to be subsampled.
273 """
275 def wrapper(x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]:
276 return transform.apply_forward(x=x, y=y)
278 return Region(
279 # Shapely overloads don't seem to have been annotated rigorously
280 shapely.transform(self._impl, wrapper, interleaved=False) # type: ignore[arg-type]
281 ).try_to_polygon()
283 def serialize(self) -> RegionSerializationModel:
284 """Serialize the region to a Pydantic model.
286 Region serialization uses a subset of the GeoJSON specification (IETF
287 RFC 7946).
288 """
289 return RegionSerializationModel.model_validate_json(shapely.to_geojson(self._impl))
291 @classmethod
292 def __get_pydantic_core_schema__(
293 cls, source_type: Any, handler: pydantic.GetCoreSchemaHandler
294 ) -> pcs.CoreSchema:
295 from_model_schema = pcs.chain_schema(
296 [
297 handler(RegionSerializationModel),
298 pcs.no_info_plain_validator_function(RegionSerializationModel.deserialize),
299 ]
300 )
301 return pcs.json_or_python_schema(
302 json_schema=from_model_schema,
303 python_schema=pcs.union_schema([pcs.is_instance_schema(cls), from_model_schema]),
304 serialization=pcs.plain_serializer_function_ser_schema(cls.serialize),
305 )
307 @classmethod
308 def __get_pydantic_json_schema__(
309 cls, schema: pcs.CoreSchema, handler: pydantic.GetJsonSchemaHandler
310 ) -> JsonSchemaValue:
311 return handler(RegionSerializationModel.__pydantic_core_schema__)
314class Polygon(Region):
315 """A simple 2-d polygon in Euclidean coordinates, with no holes.
317 Parameters
318 ----------
319 x_vertices
320 The x coordinates of the vertices of the polygon.
321 y_vertices
322 The y coordinate of the vertices of the polygon.
323 """
325 def __init__(self, *, x_vertices: npt.ArrayLike, y_vertices: npt.ArrayLike) -> None:
326 self._vertices = np.stack(
327 [np.asarray(x_vertices).flat, np.asarray(y_vertices).flat], dtype=np.float64
328 ).transpose()
329 self._vertices.flags.writeable = False
330 super().__init__(shapely.Polygon(self._vertices))
332 @staticmethod
333 def from_box(box: Box) -> Polygon:
334 """Construct from an integer-coordinate box.
336 Parameters
337 ----------
338 box
339 Integer-coordinate box to convert to a polygon.
341 Notes
342 -----
343 Because the integer min and max coordinates of the box are
344 interpreted as pixel centers, these are expanded by 0.5 on all sides
345 before using them to form the polygon vertices.
346 """
347 return Polygon(
348 x_vertices=[box.x.min - 0.5, box.x.min - 0.5, box.x.max + 0.5, box.x.max + 0.5],
349 y_vertices=[box.y.min - 0.5, box.y.max + 0.5, box.y.max + 0.5, box.y.min - 0.5],
350 )
352 @property
353 def n_vertices(self) -> int:
354 """The number of vertices in the polygon."""
355 return self._vertices.shape[0]
357 @property
358 def x_vertices(self) -> np.ndarray:
359 """The x coordinates of the vertices of the polygon.
361 This is a read-only array; polygons are immutable.
362 """
363 return self._vertices[:, 0]
365 @property
366 def y_vertices(self) -> np.ndarray:
367 """The y coordinates of the vertices of the polygon.
369 This is a read-only array; polygons are immutable.
370 """
371 return self._vertices[:, 1]
373 @property
374 def centroid(self) -> XY[float]:
375 """The centroid of the polygon (`XY` [`float`])."""
376 c = self._impl.centroid
377 return XY(x=c.x, y=c.y)
379 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report:
380 """Return a `Report` describing this polygon.
382 Parameters
383 ----------
384 options : `DescribeOptions`, optional
385 Unused; accepted for interface compatibility.
386 """
387 return Report(
388 type_name="Polygon",
389 summary=self._impl.wkt,
390 fields=[
391 ReportField(label="x_vertices", value=self.x_vertices, repr_value=repr(self.x_vertices)),
392 ReportField(label="y_vertices", value=self.y_vertices, repr_value=repr(self.y_vertices)),
393 ],
394 )
396 def __repr__(self) -> str:
397 return f"Polygon(x_vertices={self.x_vertices!r}, y_vertices={self.y_vertices!r})"
399 def transform(self, transform: Transform[Any, Any]) -> Polygon:
400 # Docstring inherited.
401 result = super().transform(transform)
402 assert isinstance(result, Polygon), "Transforming a polygon should not change its topology."
403 return result
405 @staticmethod
406 def from_legacy(legacy: LegacyPolygon) -> Polygon:
407 """Convert from a legacy `lsst.afw.geom.Polygon` instance.
409 Parameters
410 ----------
411 legacy
412 Legacy `lsst.afw.geom.Polygon` to convert.
413 """
414 vertices = legacy.getVertices()
415 x_vertices = np.zeros(len(vertices), dtype=np.float64)
416 y_vertices = np.zeros(len(vertices), dtype=np.float64)
417 for n, point in enumerate(vertices):
418 x_vertices[n] = point.x
419 y_vertices[n] = point.y
420 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices)
422 def to_legacy(self) -> LegacyPolygon:
423 """Convert to a legacy `lsst.afw.geom.Polygon` instance."""
424 from lsst.afw.geom import Polygon as LegacyPolygon
425 from lsst.geom import Point2D
427 return LegacyPolygon([Point2D(x, y) for x, y in zip(self.x_vertices, self.y_vertices)])
430class RegionSerializationModel(pydantic.BaseModel):
431 """Serialization model for `Region` and `Polygon`.
433 This model is a subset of the GeoJSON specification (IETF RFC 7946).
434 """
436 type: Literal["Polygon", "MultiPolygon"] = pydantic.Field(description="Geometry type.")
438 coordinates: list[list[tuple[float, float] | list[tuple[float, float]]]] = pydantic.Field(
439 description="Vertices of the polygon or polygons."
440 )
442 def deserialize(self) -> Region:
443 """Deserialize into a `Region` (a `Polygon`, if possible)."""
444 region_impl = shapely.from_geojson(self.model_dump_json())
445 assert isinstance(region_impl, shapely.Polygon | shapely.MultiPolygon), (
446 "Other geometry types are not used."
447 )
448 return Region(region_impl).try_to_polygon()
451if TYPE_CHECKING:
453 def _test_types() -> None:
454 region = Region(shapely.Point(0, 0).buffer(1.0))
455 arr = np.zeros(3)
457 # Region satisfies the Bounds Protocol.
458 bounds: Bounds = region
460 # Region.contains: XY/YX, scalar, array-like
461 assert_type(region.contains(x=1.0, y=2.0), bool)
462 assert_type(region.contains(x=arr, y=arr), np.ndarray)
463 assert_type(region.contains(XY(1.0, 2.0)), bool)
464 assert_type(region.contains(YX(2.0, 1.0)), bool)
465 assert_type(region.contains(XY(arr, arr)), np.ndarray)
466 assert_type(region.contains(YX(arr, arr)), np.ndarray)
468 # Via the Bounds Protocol view, same signatures hold.
469 assert_type(bounds.contains(x=1, y=1), bool)
470 assert_type(bounds.contains(x=arr, y=arr), np.ndarray)
471 assert_type(bounds.contains(XY(1, 1)), bool)
472 assert_type(bounds.contains(YX(1, 1)), bool)
473 assert_type(bounds.contains(XY(arr, arr)), np.ndarray)
474 assert_type(bounds.contains(YX(arr, arr)), np.ndarray)