Coverage for tests/test_polygon.py: 99%

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 

14import dataclasses 

15 

16import astropy.units as u 

17import numpy as np 

18import pydantic 

19import pytest 

20import shapely 

21 

22from lsst.images import ( 

23 XY, 

24 YX, 

25 Box, 

26 GeneralFrame, 

27 NoOverlapError, 

28 Polygon, 

29 Region, 

30 RegionDomainError, 

31 RegionSerializationModel, 

32 Transform, 

33) 

34from lsst.images.tests import assert_values_equal, check_bounds_contains_broadcasting 

35 

36try: 

37 import lsst.afw.geom # noqa: F401 

38 

39 have_legacy = True 

40except ImportError: 

41 have_legacy = False 

42 

43skip_no_legacy = pytest.mark.skipif(not have_legacy, reason="lsst legacy packages could not be imported.") 

44 

45 

46class _PolygonHolder(pydantic.BaseModel): 

47 polygon: Polygon 

48 

49 

50class _RegionHolder(pydantic.BaseModel): 

51 region: Region 

52 

53 

54def _make_polygon() -> Polygon: 

55 """Return a near-box quadrilateral that is easy to reason about.""" 

56 x_vertices = [32.0, 31.0, 50.0, 53.0] 

57 y_vertices = [-5.0, 7.0, 7.2, -4.8] 

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

59 

60 

61@dataclasses.dataclass 

62class _TestRegions: 

63 """Four overlapping box-polygons for region operation tests. 

64 

65 Rough layout (y increasing upwards):: 

66 

67 ┌─────┐ 

68 ┌───┼─┐B┌─┼────┐ 

69 │ A└─┼─┼─┘ ┌─┐│ 

70 └─────┘ │ C │D││ 

71 │ └─┘│ 

72 └──────┘ 

73 """ 

74 

75 a: Polygon 

76 b: Polygon 

77 c: Polygon 

78 d: Polygon 

79 

80 def __init__(self) -> None: 

81 self.a = Polygon.from_box(Box.factory[3:6, 0:5]) 

82 self.b = Polygon.from_box(Box.factory[4:7, 3:8]) 

83 self.c = Polygon.from_box(Box.factory[0:6, 6:12]) 

84 self.d = Polygon.from_box(Box.factory[1:4, 9:10]) 

85 

86 

87def test_vertices() -> None: 

88 """Test the vertices accessors.""" 

89 polygon = _make_polygon() 

90 x_vertices = [32.0, 31.0, 50.0, 53.0] 

91 y_vertices = [-5.0, 7.0, 7.2, -4.8] 

92 assert polygon.n_vertices == 4 

93 np.testing.assert_array_equal(polygon.x_vertices, np.asarray(x_vertices)) 

94 np.testing.assert_array_equal(polygon.y_vertices, np.asarray(y_vertices)) 

95 with pytest.raises(ValueError, match="read-only"): 

96 polygon.x_vertices[0] = 0.0 

97 with pytest.raises(ValueError, match="read-only"): 

98 polygon.y_vertices[0] = 0.0 

99 

100 

101def test_boxes() -> None: 

102 """Test 'from_box', the ``area`` property, and the 'contains' method 

103 with polygon arguments. 

104 """ 

105 polygon = _make_polygon() 

106 small = Polygon.from_box(Box.factory[-3:3, 40:45]) 

107 assert small.area == 30.0 

108 assert small.bbox == Box.factory[-3:3, 40:45] 

109 assert polygon.contains(small) 

110 assert not small.contains(polygon) 

111 assert small.centroid == XY(x=42.0, y=-0.5) 

112 big = Polygon.from_box(Box.factory[-10:10, 20:60]) 

113 assert big.area == 800.0 

114 assert not polygon.contains(big) 

115 assert big.contains(polygon) 

116 medium = Polygon.from_box(Box.factory[-4:8, 31:52]) 

117 assert medium.area == 252.0 

118 assert not polygon.contains(medium) 

119 assert not medium.contains(polygon) 

120 assert polygon.contains(polygon) 

121 

122 

123def test_transform() -> None: 

124 """Test applying a coordinate transform to a polygon.""" 

125 polygon = _make_polygon() 

126 matrix = np.array([[0.4, 0.25], [-0.20, 0.6]]) 

127 t = Transform.affine(GeneralFrame(unit=u.pix), GeneralFrame(unit=u.pix), matrix) 

128 tp = polygon.transform(t) 

129 assert isinstance(tp, Polygon) 

130 assert_values_equal(tp.area, polygon.area * np.linalg.det(matrix), rtol=1e-5) 

131 xyt = t.apply_forward(x=polygon.x_vertices, y=polygon.y_vertices) 

132 # Slicing below is because shapely sometimes adds a duplicate closing 

133 # vertex. 

134 assert_values_equal(tp.x_vertices[: len(xyt.x)], xyt.x, rtol=1e-14) 

135 assert_values_equal(tp.y_vertices[: len(xyt.y)], xyt.y, rtol=1e-14) 

136 

137 

138def test_transform_outside_domain() -> None: 

139 """Test that a transform yielding non-finite vertices raises a 

140 `RegionDomainError` instead of producing a broken region. 

141 """ 

142 polygon = _make_polygon() 

143 matrix = np.array([[1e308, 0.0], [0.0, 1e308]]) 

144 t = Transform.affine(GeneralFrame(unit=u.pix), GeneralFrame(unit=u.pix), matrix) 

145 assert not np.isfinite(t.apply_forward(x=polygon.x_vertices, y=polygon.y_vertices).x).all() 

146 with pytest.raises(RegionDomainError): 

147 polygon.transform(t) 

148 

149 

150def test_contains_points() -> None: 

151 """Test the 'contains' method with points.""" 

152 polygon = _make_polygon() 

153 assert polygon.contains(x=40.0, y=0.0) 

154 assert not polygon.contains(x=0.0, y=0.0) 

155 assert not polygon.contains(x=40.0, y=10.0) 

156 np.testing.assert_array_equal( 

157 polygon.contains(x=np.array([40.0, 0.0, 40.0]), y=np.array([0.0, 0.0, 10.0])), 

158 np.array([True, False, False]), 

159 ) 

160 

161 

162def test_contains_points_xy_yx() -> None: 

163 """Verify that Region.contains accepts XY and YX positional arguments.""" 

164 polygon = _make_polygon() 

165 # Scalar XY and YX — results must match the keyword form. 

166 assert polygon.contains(XY(x=40.0, y=0.0)) == polygon.contains(x=40.0, y=0.0) 

167 assert polygon.contains(YX(y=0.0, x=40.0)) == polygon.contains(x=40.0, y=0.0) 

168 assert not polygon.contains(XY(x=0.0, y=0.0)) 

169 assert not polygon.contains(YX(y=0.0, x=0.0)) 

170 # Array XY and YX. 

171 xv = np.array([40.0, 0.0, 40.0]) 

172 yv = np.array([0.0, 0.0, 10.0]) 

173 np.testing.assert_array_equal(polygon.contains(XY(xv, yv)), polygon.contains(x=xv, y=yv)) 

174 np.testing.assert_array_equal(polygon.contains(YX(yv, xv)), polygon.contains(x=xv, y=yv)) 

175 # Mixing positional point with keyword x= or y= must raise TypeError. 

176 with pytest.raises(TypeError): 

177 polygon.contains(XY(x=40.0, y=0.0), x=40.0) 

178 with pytest.raises(TypeError): 

179 polygon.contains(YX(y=0.0, x=40.0), y=0.0) 

180 

181 

182def test_region_contains_broadcasting() -> None: 

183 """Test that Region.contains broadcasts like a numpy ufunc.""" 

184 check_bounds_contains_broadcasting(_make_polygon()) 

185 

186 

187def test_polygon_io() -> None: 

188 """Test serialization and stringification.""" 

189 polygon = _make_polygon() 

190 assert ( 

191 RegionSerializationModel.model_validate_json(polygon.serialize().model_dump_json()).deserialize() 

192 == polygon 

193 ) 

194 assert Polygon.from_wkt(polygon.wkt) == polygon 

195 assert Polygon.from_wkt(str(polygon)) == polygon 

196 assert eval(repr(polygon), {"array": np.array, "Polygon": Polygon}) == polygon 

197 

198 

199def test_polygon_model_field() -> None: 

200 """Test that we can use a Polygon directly as a Pydantic model field.""" 

201 polygon = _make_polygon() 

202 holder = _PolygonHolder(polygon=polygon) 

203 assert polygon == holder.model_validate_json(holder.model_dump_json()).polygon 

204 assert ( 

205 _PolygonHolder.model_json_schema()["properties"]["polygon"] 

206 == RegionSerializationModel.model_json_schema() 

207 ) 

208 

209 

210@skip_no_legacy 

211def test_polygon_legacy() -> None: 

212 """Test conversion to/from lsst.afw.geom.Polygon.""" 

213 polygon = _make_polygon() 

214 legacy_polygon = polygon.to_legacy() 

215 assert legacy_polygon.calculateArea() == polygon.area 

216 assert Polygon.from_legacy(legacy_polygon) == polygon 

217 

218 

219def test_intersection() -> None: 

220 """Test region intersection.""" 

221 regions = _TestRegions() 

222 # Usual case: 

223 assert regions.a.intersection(regions.b) == Polygon.from_box(Box.factory[4:6, 3:5]) 

224 # No-overlap case: 

225 with pytest.raises(NoOverlapError): 

226 regions.a.intersection(regions.c) 

227 # LHS fully contains RHS: 

228 assert regions.c.intersection(regions.d) == regions.d 

229 # Intersections with the boxes themselves should return boxes: 

230 assert regions.a.intersection(regions.b.bbox) == Box.factory[4:6, 3:5] 

231 assert regions.a.bbox.intersection(regions.b) == Box.factory[4:6, 3:5] 

232 assert ( 

233 # A Box is not possible when the result is not simple. 

234 regions.a.union(regions.c).intersection(regions.b.bbox) 

235 == regions.a.union(regions.c).intersection(regions.b) 

236 ) 

237 

238 

239def test_union() -> None: 

240 """Test region union.""" 

241 regions = _TestRegions() 

242 # Usual case: 

243 assert regions.a.union(regions.b).bbox == Box.factory[3:7, 0:8] 

244 assert regions.a.union(regions.b).area == 15 + 15 - 4 

245 # Operands are disjoint, so union is not a single Polygon: 

246 assert not isinstance(regions.a.union(regions.c), Polygon) 

247 assert regions.a.union(regions.c).area == regions.a.area + regions.c.area 

248 # LHS fully contains RHS: 

249 assert regions.c.union(regions.d) == regions.c 

250 

251 

252def test_difference() -> None: 

253 """Test region difference.""" 

254 regions = _TestRegions() 

255 # Usual case: 

256 assert regions.a.difference(regions.b).bbox == regions.a.bbox 

257 assert regions.a.difference(regions.b).area == 15 - 4 

258 # Operands are disjoint, so difference is just the LHS. 

259 assert regions.a.difference(regions.c) == regions.a 

260 # LHS fully contains RHS -> polygon with hole is not a Polygon: 

261 assert not isinstance(regions.c.difference(regions.d), Polygon) 

262 assert regions.c.difference(regions.d).bbox == regions.c.bbox 

263 assert regions.c.difference(regions.d).area == regions.c.area - regions.d.area 

264 # RHS fully contains LHS -> region is empty. 

265 assert regions.d.difference(regions.d).area == 0 

266 

267 

268def test_region_io() -> None: 

269 """Test serialization and stringification of non-polygon regions.""" 

270 regions = _TestRegions() 

271 # A two-polygon region with a hole: 

272 region = regions.a.union(regions.c).difference(regions.d) 

273 assert ( 

274 RegionSerializationModel.model_validate_json(region.serialize().model_dump_json()).deserialize() 

275 == region 

276 ) 

277 assert Region.from_wkt(region.wkt) == region 

278 assert Region.from_wkt(str(region)) == region 

279 assert eval(repr(region), {"Region": Region}) == region 

280 

281 

282def test_region_model_field() -> None: 

283 """Test that we can use a Region directly as a Pydantic model field.""" 

284 regions = _TestRegions() 

285 region = regions.a.union(regions.c).difference(regions.d) 

286 holder = _RegionHolder(region=region) 

287 assert region == holder.model_validate_json(holder.model_dump_json()).region 

288 assert ( 

289 _RegionHolder.model_json_schema()["properties"]["region"] 

290 == RegionSerializationModel.model_json_schema() 

291 ) 

292 

293 

294def test_region_polygon_repr_str_pinned() -> None: 

295 """Region and Polygon str/repr match their documented forms.""" 

296 poly = Polygon(x_vertices=[0.0, 4.0, 4.0, 0.0], y_vertices=[0.0, 0.0, 5.0, 5.0]) 

297 assert repr(poly) == f"Polygon(x_vertices={poly.x_vertices!r}, y_vertices={poly.y_vertices!r})" 

298 # str is inherited from Region and uses WKT. 

299 assert str(poly) == poly.wkt 

300 

301 # Region.from_wkt on a simple polygon returns a Polygon via try_to_polygon. 

302 # To test Region.__str__/__repr__ directly, construct a true Region 

303 # (multi-polygon), which is not coerced to Polygon. 

304 multi = shapely.MultiPolygon( 

305 [ 

306 shapely.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]), 

307 shapely.Polygon([(2, 2), (3, 2), (3, 3), (2, 3)]), 

308 ] 

309 ) 

310 region = Region(multi) 

311 assert str(region) == region.wkt 

312 assert repr(region) == f"Region.from_wkt({region.wkt!r})"