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 RegionSerializationModel, 

31 Transform, 

32) 

33from lsst.images.tests import assert_close, check_bounds_contains_broadcasting 

34 

35try: 

36 import lsst.afw.geom # noqa: F401 

37 

38 have_legacy = True 

39except ImportError: 

40 have_legacy = False 

41 

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

43 

44 

45class _PolygonHolder(pydantic.BaseModel): 

46 polygon: Polygon 

47 

48 

49class _RegionHolder(pydantic.BaseModel): 

50 region: Region 

51 

52 

53def _make_polygon() -> Polygon: 

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

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

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

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

58 

59 

60@dataclasses.dataclass 

61class _TestRegions: 

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

63 

64 Rough layout (y increasing upwards):: 

65 

66 ┌─────┐ 

67 ┌───┼─┐B┌─┼────┐ 

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

69 └─────┘ │ C │D││ 

70 │ └─┘│ 

71 └──────┘ 

72 """ 

73 

74 a: Polygon 

75 b: Polygon 

76 c: Polygon 

77 d: Polygon 

78 

79 def __init__(self) -> None: 

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

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

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

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

84 

85 

86def test_vertices() -> None: 

87 """Test the vertices accessors.""" 

88 polygon = _make_polygon() 

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

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

91 assert polygon.n_vertices == 4 

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

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

94 with pytest.raises(ValueError): 

95 polygon.x_vertices[0] = 0.0 

96 with pytest.raises(ValueError): 

97 polygon.y_vertices[0] = 0.0 

98 

99 

100def test_boxes() -> None: 

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

102 with polygon arguments. 

103 """ 

104 polygon = _make_polygon() 

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

106 assert small.area == 30.0 

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

108 assert polygon.contains(small) 

109 assert not small.contains(polygon) 

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

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

112 assert big.area == 800.0 

113 assert not polygon.contains(big) 

114 assert big.contains(polygon) 

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

116 assert medium.area == 252.0 

117 assert not polygon.contains(medium) 

118 assert not medium.contains(polygon) 

119 assert polygon.contains(polygon) 

120 

121 

122def test_transform() -> None: 

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

124 polygon = _make_polygon() 

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

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

127 tp = polygon.transform(t) 

128 assert isinstance(tp, Polygon) 

129 assert_close(tp.area, polygon.area * np.linalg.det(matrix)) 

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

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

132 # vertex. 

133 assert_close(tp.x_vertices[: len(xyt.x)], xyt.x) 

134 assert_close(tp.y_vertices[: len(xyt.y)], xyt.y) 

135 

136 

137def test_contains_points() -> None: 

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

139 polygon = _make_polygon() 

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

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

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

143 np.testing.assert_array_equal( 

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

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

146 ) 

147 

148 

149def test_contains_points_xy_yx() -> None: 

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

151 polygon = _make_polygon() 

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

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

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

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

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

157 # Array XY and YX. 

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

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

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

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

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

163 with pytest.raises(TypeError): 

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

165 with pytest.raises(TypeError): 

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

167 

168 

169def test_region_contains_broadcasting() -> None: 

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

171 check_bounds_contains_broadcasting(_make_polygon()) 

172 

173 

174def test_polygon_io() -> None: 

175 """Test serialization and stringification.""" 

176 polygon = _make_polygon() 

177 assert ( 

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

179 == polygon 

180 ) 

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

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

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

184 

185 

186def test_polygon_model_field() -> None: 

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

188 polygon = _make_polygon() 

189 holder = _PolygonHolder(polygon=polygon) 

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

191 assert ( 

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

193 == RegionSerializationModel.model_json_schema() 

194 ) 

195 

196 

197@skip_no_legacy 

198def test_polygon_legacy() -> None: 

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

200 polygon = _make_polygon() 

201 legacy_polygon = polygon.to_legacy() 

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

203 assert Polygon.from_legacy(legacy_polygon) == polygon 

204 

205 

206def test_intersection() -> None: 

207 """Test region intersection.""" 

208 regions = _TestRegions() 

209 # Usual case: 

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

211 # No-overlap case: 

212 with pytest.raises(NoOverlapError): 

213 regions.a.intersection(regions.c) 

214 # LHS fully contains RHS: 

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

216 # Intersections with the boxes themselves should return boxes: 

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

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

219 assert ( 

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

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

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

223 ) 

224 

225 

226def test_union() -> None: 

227 """Test region union.""" 

228 regions = _TestRegions() 

229 # Usual case: 

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

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

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

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

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

235 # LHS fully contains RHS: 

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

237 

238 

239def test_difference() -> None: 

240 """Test region difference.""" 

241 regions = _TestRegions() 

242 # Usual case: 

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

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

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

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

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

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

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

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

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

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

253 

254 

255def test_region_io() -> None: 

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

257 regions = _TestRegions() 

258 # A two-polygon region with a hole: 

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

260 assert ( 

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

262 == region 

263 ) 

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

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

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

267 

268 

269def test_region_model_field() -> None: 

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

271 regions = _TestRegions() 

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

273 holder = _RegionHolder(region=region) 

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

275 assert ( 

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

277 == RegionSerializationModel.model_json_schema() 

278 ) 

279 

280 

281def test_region_polygon_repr_str_pinned() -> None: 

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

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

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

285 # str is inherited from Region and uses WKT. 

286 assert str(poly) == poly.wkt 

287 

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

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

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

291 multi = shapely.MultiPolygon( 

292 [ 

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

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

295 ] 

296 ) 

297 region = Region(multi) 

298 assert str(region) == region.wkt 

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