Coverage for python/lsst/images/_polygon.py: 49%

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. 

11 

12from __future__ import annotations 

13 

14__all__ = ("Polygon", "Region", "RegionDomainError", "RegionSerializationModel") 

15 

16from typing import TYPE_CHECKING, Any, Literal, assert_type, overload 

17 

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 

24 

25from ._geom import XY, YX, Bounds, Box 

26from .describe import DescribableMixin, DescribeOptions, Report, ReportField 

27from .utils import round_half_down, round_half_up 

28 

29if TYPE_CHECKING: 

30 from ._geom import Bounds 

31 from ._transforms import Transform 

32 

33 try: 

34 from lsst.afw.geom import Polygon as LegacyPolygon 

35 except ImportError: 

36 type LegacyPolygon = Any # type: ignore[no-redef] 

37 

38 

39class RegionDomainError(ValueError): 

40 """Exception raised when a region is transformed to non-finite 

41 coordinates, i.e. it lies outside the domain of the transform. 

42 """ 

43 

44 

45class Region(DescribableMixin): 

46 """A 2-d Euclidean region represented as one or more polygons with 

47 optional holes. 

48 

49 Parameters 

50 ---------- 

51 geometry 

52 A polygon or multi-polygon from the Shapely library. 

53 """ 

54 

55 def __init__(self, geometry: shapely.Polygon | shapely.MultiPolygon) -> None: 

56 self._impl = geometry 

57 

58 @property 

59 def area(self) -> float: 

60 """The area of the region (`float`).""" 

61 return self._impl.area 

62 

63 @property 

64 def bbox(self) -> Box: 

65 """The integer-coordinate bounding box of the region (`Box`). 

66 

67 Because a `Box` logically contains the entirety of the pixels on its 

68 boundary, but the centers of those pixels are the numerical values of 

69 its bounds, the region may contain vertices that are up to 0.5 beyond 

70 the integer box coordinates in either dimension. 

71 """ 

72 x_min, y_min, x_max, y_max = self._impl.bounds 

73 return Box.factory[ 

74 round_half_up(y_min) : round_half_down(y_max) + 1, 

75 round_half_up(x_min) : round_half_down(x_max) + 1, 

76 ] 

77 

78 @property 

79 def wkt(self) -> str: 

80 """The 'Well-Known Text' representation of this region (`str`).""" 

81 return self._impl.wkt 

82 

83 @staticmethod 

84 def from_wkt(wkt: str) -> Region: 

85 """Construct from a 'Well-Known Text' string. 

86 

87 Parameters 

88 ---------- 

89 wkt 

90 Well-Known Text representation of the region. 

91 """ 

92 impl = shapely.from_wkt(wkt) 

93 if not isinstance(impl, shapely.Polygon | shapely.MultiPolygon): 93 ↛ 94line 93 didn't jump to line 94 because the condition on line 93 was never true

94 raise ValueError("Only Polygon and MultiPolygon geometries can be converted to Regions.") 

95 return Region(impl).try_to_polygon() 

96 

97 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report: 

98 """Return a `Report` describing this region. 

99 

100 Parameters 

101 ---------- 

102 options : `DescribeOptions`, optional 

103 Unused; accepted for interface compatibility. 

104 """ 

105 return Report( 

106 type_name="Region", 

107 summary=self._impl.wkt, 

108 fields=[ 

109 ReportField( 

110 label="wkt", 

111 value=self._impl.wkt, 

112 repr_value=repr(self._impl.wkt), 

113 positional=True, 

114 ) 

115 ], 

116 ) 

117 

118 def __repr__(self) -> str: 

119 return f"Region.from_wkt({self._impl.wkt!r})" 

120 

121 def __eq__(self, other: object) -> bool: 

122 if isinstance(other, Region): 

123 return bool(shapely.equals(self._impl, other._impl)) 

124 return NotImplemented 

125 

126 @overload 

127 def contains(self, point: XY[int | float] | YX[int | float], /) -> bool: ... 

128 

129 @overload 

130 def contains(self, point: XY[npt.ArrayLike] | YX[npt.ArrayLike], /) -> np.ndarray: ... 

131 

132 @overload 

133 def contains(self, /, *, x: int | float, y: int | float) -> bool: ... 

134 

135 @overload 

136 def contains(self, /, *, x: npt.ArrayLike, y: npt.ArrayLike) -> np.ndarray: ... 

137 

138 @overload 

139 def contains(self, other: Polygon, /) -> bool: ... 

140 

141 def contains( 

142 self, 

143 other: Region | XY[Any] | YX[Any] | None = None, 

144 /, 

145 *, 

146 x: Any = None, 

147 y: Any = None, 

148 ) -> bool | np.ndarray: 

149 """Test whether the geometry contains the given points or another 

150 geometry. 

151 

152 Parameters 

153 ---------- 

154 other 

155 Another geometry, or an `XY` or `YX` coordinate pair, to compare 

156 to. Mutually exclusive with ``x`` and ``y``. 

157 x 

158 One or more X coordinates to test for containment, as a scalar or 

159 any array-like. Results are broadcast against ``y``. 

160 Mutually exclusive with ``other``. 

161 y 

162 One or more Y coordinates to test for containment, as a scalar or 

163 any array-like. Results are broadcast against ``x``. 

164 Mutually exclusive with ``other``. 

165 """ 

166 match other: 

167 case None: 

168 if x is None or y is None: 168 ↛ 169line 168 didn't jump to line 169 because the condition on line 168 was never true

169 raise TypeError("Pass either a point or both x= and y= to 'Region.contains'.") 

170 case Region(): 

171 if x is not None or y is not None: 171 ↛ 172line 171 didn't jump to line 172 because the condition on line 171 was never true

172 raise TypeError("'Region.contains' point argument is mutually exclusive with x= and y=.") 

173 return self._impl.contains(other._impl) 

174 case XY() | YX(): 174 ↛ 178line 174 didn't jump to line 178 because the pattern on line 174 always matched

175 if x is not None or y is not None: 

176 raise TypeError("'Region.contains' point argument is mutually exclusive with x= and y=.") 

177 x, y = other.x, other.y 

178 case _: 

179 raise TypeError(f"Unexpected positional argument type: {type(other)!r}.") 

180 return shapely.contains_xy(self._impl, x=x, y=y) 

181 

182 @overload 

183 def intersection(self, other: Region) -> Region: ... 

184 

185 @overload 

186 def intersection(self, other: Box) -> Region | Box: ... 

187 

188 @overload 

189 def intersection(self, other: Bounds) -> Bounds: ... 

190 

191 def intersection(self, other: Bounds) -> Bounds: 

192 """Compute the intersection of this region with a `Bounds` object. 

193 

194 Notes 

195 ----- 

196 Because `Region` implements the `Bounds` interface, its intersections 

197 need to support all other `Bounds` objects. This is not true of other 

198 `Region` point-set operations like `union` and `difference`. 

199 

200 Parameters 

201 ---------- 

202 other 

203 Bounds to intersect with this region. 

204 """ 

205 from ._concrete_bounds import _intersect_region 

206 

207 return _intersect_region(self, other) 

208 

209 def union(self, other: Region) -> Region: 

210 """Compute the point-set union of this region with another. 

211 

212 Parameters 

213 ---------- 

214 other 

215 Region to union with this one. 

216 """ 

217 impl = shapely.union(self._impl, other._impl) 

218 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), ( 

219 "A union of Polygons and MultiPolygons should be one of those." 

220 ) 

221 return Region(impl).try_to_polygon() 

222 

223 def difference(self, other: Region) -> Region: 

224 """Compute the point-set difference of this region with another. 

225 

226 Parameters 

227 ---------- 

228 other 

229 Region to subtract from this one. 

230 """ 

231 impl = shapely.difference(self._impl, other._impl) 

232 assert isinstance(impl, shapely.Polygon | shapely.MultiPolygon), ( 

233 "A difference of Polygons and MultiPolygons should be one of those." 

234 ) 

235 return Region(impl).try_to_polygon() 

236 

237 def try_to_polygon(self) -> Region: 

238 """If the underlying geometry is a single polygon with no holes, 

239 return a `Polygon` instance holding it. 

240 

241 In all other cases ``self`` is returned. 

242 """ 

243 impl = self._impl 

244 if isinstance(impl, shapely.MultiPolygon) and len(impl.geoms) == 1: 244 ↛ 245line 244 didn't jump to line 245 because the condition on line 244 was never true

245 impl = impl.geoms[0] 

246 if isinstance(impl, shapely.Polygon) and not impl.interiors: 

247 vertices = np.array(impl.exterior.coords) 

248 return Polygon(x_vertices=vertices[:, 0], y_vertices=vertices[:, 1]) 

249 return self 

250 

251 def try_to_box(self) -> Region | Box: 

252 """If the underlying geometry is a rectangle that fully covers integer 

253 pixels (i.e. has all vertices at half-integer positions), return the 

254 equivalent `Box`. 

255 

256 In all other cases ``self`` is returned. 

257 """ 

258 if Polygon.from_box(self.bbox) == self: 

259 return self.bbox 

260 return self 

261 

262 def to_shapely(self) -> shapely.Polygon | shapely.MultiPolygon: 

263 """Convert to a `shapely.Polygon` or `shapely.MultiPolygon` object.""" 

264 return self._impl 

265 

266 def transform(self, transform: Transform[Any, Any]) -> Region: 

267 """Transform the coordinates of the region, returning a new one. 

268 

269 Parameters 

270 ---------- 

271 transform 

272 Coordinate transform to apply (in the forward direction). 

273 

274 Returns 

275 ------- 

276 `Region` 

277 The transformed region, as a `Polygon` if it is a single polygon 

278 with no holes. 

279 

280 Raises 

281 ------ 

282 RegionDomainError 

283 Raised if any transformed vertex is not finite, as happens when 

284 the region lies outside the domain of a projection. 

285 

286 Notes 

287 ----- 

288 This applies the transform to all vertices, assuming that the 

289 transform is close enough to affine that the topology of the geometry 

290 does not change and straight-line edges do not need to be subsampled. 

291 """ 

292 

293 def wrapper(x: np.ndarray, y: np.ndarray) -> XY[np.ndarray]: 

294 result = transform.apply_forward(x=x, y=y) 

295 if not (np.isfinite(result.x).all() and np.isfinite(result.y).all()): 

296 raise RegionDomainError(f"Region {self} transforms to non-finite coordinates.") 

297 return result 

298 

299 return Region( 

300 # Shapely overloads don't seem to have been annotated rigorously 

301 shapely.transform(self._impl, wrapper, interleaved=False) # type: ignore[arg-type] 

302 ).try_to_polygon() 

303 

304 def serialize(self) -> RegionSerializationModel: 

305 """Serialize the region to a Pydantic model. 

306 

307 Region serialization uses a subset of the GeoJSON specification (IETF 

308 RFC 7946). 

309 """ 

310 return RegionSerializationModel.model_validate_json(shapely.to_geojson(self._impl)) 

311 

312 @classmethod 

313 def __get_pydantic_core_schema__( 

314 cls, source_type: Any, handler: pydantic.GetCoreSchemaHandler 

315 ) -> pcs.CoreSchema: 

316 from_model_schema = pcs.chain_schema( 

317 [ 

318 handler(RegionSerializationModel), 

319 pcs.no_info_plain_validator_function(RegionSerializationModel.deserialize), 

320 ] 

321 ) 

322 return pcs.json_or_python_schema( 

323 json_schema=from_model_schema, 

324 python_schema=pcs.union_schema([pcs.is_instance_schema(cls), from_model_schema]), 

325 serialization=pcs.plain_serializer_function_ser_schema(cls.serialize), 

326 ) 

327 

328 @classmethod 

329 def __get_pydantic_json_schema__( 

330 cls, schema: pcs.CoreSchema, handler: pydantic.GetJsonSchemaHandler 

331 ) -> JsonSchemaValue: 

332 return handler(RegionSerializationModel.__pydantic_core_schema__) 

333 

334 

335class Polygon(Region): 

336 """A simple 2-d polygon in Euclidean coordinates, with no holes. 

337 

338 Parameters 

339 ---------- 

340 x_vertices 

341 The x coordinates of the vertices of the polygon. 

342 y_vertices 

343 The y coordinate of the vertices of the polygon. 

344 """ 

345 

346 def __init__(self, *, x_vertices: npt.ArrayLike, y_vertices: npt.ArrayLike) -> None: 

347 self._vertices = np.stack( 

348 [np.asarray(x_vertices).flat, np.asarray(y_vertices).flat], dtype=np.float64 

349 ).transpose() 

350 self._vertices.flags.writeable = False 

351 super().__init__(shapely.Polygon(self._vertices)) 

352 

353 @staticmethod 

354 def from_box(box: Box) -> Polygon: 

355 """Construct from an integer-coordinate box. 

356 

357 Parameters 

358 ---------- 

359 box 

360 Integer-coordinate box to convert to a polygon. 

361 

362 Notes 

363 ----- 

364 Because the integer min and max coordinates of the box are 

365 interpreted as pixel centers, these are expanded by 0.5 on all sides 

366 before using them to form the polygon vertices. 

367 """ 

368 return Polygon( 

369 x_vertices=[box.x.min - 0.5, box.x.min - 0.5, box.x.max + 0.5, box.x.max + 0.5], 

370 y_vertices=[box.y.min - 0.5, box.y.max + 0.5, box.y.max + 0.5, box.y.min - 0.5], 

371 ) 

372 

373 @property 

374 def n_vertices(self) -> int: 

375 """The number of vertices in the polygon.""" 

376 return self._vertices.shape[0] 

377 

378 @property 

379 def x_vertices(self) -> np.ndarray: 

380 """The x coordinates of the vertices of the polygon. 

381 

382 This is a read-only array; polygons are immutable. 

383 """ 

384 return self._vertices[:, 0] 

385 

386 @property 

387 def y_vertices(self) -> np.ndarray: 

388 """The y coordinates of the vertices of the polygon. 

389 

390 This is a read-only array; polygons are immutable. 

391 """ 

392 return self._vertices[:, 1] 

393 

394 @property 

395 def centroid(self) -> XY[float]: 

396 """The centroid of the polygon (`XY` [`float`]).""" 

397 c = self._impl.centroid 

398 return XY(x=c.x, y=c.y) 

399 

400 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report: 

401 """Return a `Report` describing this polygon. 

402 

403 Parameters 

404 ---------- 

405 options : `DescribeOptions`, optional 

406 Unused; accepted for interface compatibility. 

407 """ 

408 return Report( 

409 type_name="Polygon", 

410 summary=self._impl.wkt, 

411 fields=[ 

412 ReportField(label="x_vertices", value=self.x_vertices, repr_value=repr(self.x_vertices)), 

413 ReportField(label="y_vertices", value=self.y_vertices, repr_value=repr(self.y_vertices)), 

414 ], 

415 ) 

416 

417 def __repr__(self) -> str: 

418 return f"Polygon(x_vertices={self.x_vertices!r}, y_vertices={self.y_vertices!r})" 

419 

420 def transform(self, transform: Transform[Any, Any]) -> Polygon: 

421 # Docstring inherited. 

422 result = super().transform(transform) 

423 assert isinstance(result, Polygon), "Transforming a polygon should not change its topology." 

424 return result 

425 

426 @staticmethod 

427 def from_legacy(legacy: LegacyPolygon) -> Polygon: 

428 """Convert from a legacy `lsst.afw.geom.Polygon` instance. 

429 

430 Parameters 

431 ---------- 

432 legacy 

433 Legacy `lsst.afw.geom.Polygon` to convert. 

434 """ 

435 vertices = legacy.getVertices() 

436 x_vertices = np.zeros(len(vertices), dtype=np.float64) 

437 y_vertices = np.zeros(len(vertices), dtype=np.float64) 

438 for n, point in enumerate(vertices): 

439 x_vertices[n] = point.x 

440 y_vertices[n] = point.y 

441 return Polygon(x_vertices=x_vertices, y_vertices=y_vertices) 

442 

443 def to_legacy(self) -> LegacyPolygon: 

444 """Convert to a legacy `lsst.afw.geom.Polygon` instance.""" 

445 from lsst.afw.geom import Polygon as LegacyPolygon 

446 from lsst.geom import Point2D 

447 

448 return LegacyPolygon([Point2D(x, y) for x, y in zip(self.x_vertices, self.y_vertices)]) 

449 

450 

451class RegionSerializationModel(pydantic.BaseModel): 

452 """Serialization model for `Region` and `Polygon`. 

453 

454 This model is a subset of the GeoJSON specification (IETF RFC 7946). 

455 """ 

456 

457 type: Literal["Polygon", "MultiPolygon"] = pydantic.Field(description="Geometry type.") 

458 

459 coordinates: list[list[tuple[float, float] | list[tuple[float, float]]]] = pydantic.Field( 

460 description="Vertices of the polygon or polygons." 

461 ) 

462 

463 def deserialize(self) -> Region: 

464 """Deserialize into a `Region` (a `Polygon`, if possible).""" 

465 region_impl = shapely.from_geojson(self.model_dump_json()) 

466 assert isinstance(region_impl, shapely.Polygon | shapely.MultiPolygon), ( 

467 "Other geometry types are not used." 

468 ) 

469 return Region(region_impl).try_to_polygon() 

470 

471 

472if TYPE_CHECKING: 

473 

474 def _test_types() -> None: 

475 region = Region(shapely.Point(0, 0).buffer(1.0)) 

476 arr = np.zeros(3) 

477 

478 # Region satisfies the Bounds Protocol. 

479 bounds: Bounds = region 

480 

481 # Region.contains: XY/YX, scalar, array-like 

482 assert_type(region.contains(x=1.0, y=2.0), bool) 

483 assert_type(region.contains(x=arr, y=arr), np.ndarray) 

484 assert_type(region.contains(XY(1.0, 2.0)), bool) 

485 assert_type(region.contains(YX(2.0, 1.0)), bool) 

486 assert_type(region.contains(XY(arr, arr)), np.ndarray) 

487 assert_type(region.contains(YX(arr, arr)), np.ndarray) 

488 

489 # Via the Bounds Protocol view, same signatures hold. 

490 assert_type(bounds.contains(x=1, y=1), bool) 

491 assert_type(bounds.contains(x=arr, y=arr), np.ndarray) 

492 assert_type(bounds.contains(XY(1, 1)), bool) 

493 assert_type(bounds.contains(YX(1, 1)), bool) 

494 assert_type(bounds.contains(XY(arr, arr)), np.ndarray) 

495 assert_type(bounds.contains(YX(arr, arr)), np.ndarray)