test_set_operations.py 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481
  1. import numpy as np
  2. import pytest
  3. import shapely
  4. from shapely import Geometry, GeometryCollection, Polygon
  5. from shapely.testing import assert_geometries_equal
  6. from shapely.tests.common import all_types, empty, ignore_invalid, point, polygon
  7. pytestmark = pytest.mark.filterwarnings(
  8. "ignore:The symmetric_difference_all function:DeprecationWarning"
  9. )
  10. # fixed-precision operations raise GEOS exceptions on mixed dimension geometry
  11. # collections
  12. all_single_types = np.array(all_types)[
  13. ~shapely.is_empty(all_types)
  14. & (shapely.get_type_id(all_types) != shapely.GeometryType.GEOMETRYCOLLECTION)
  15. ]
  16. SET_OPERATIONS = (
  17. shapely.difference,
  18. shapely.intersection,
  19. shapely.symmetric_difference,
  20. shapely.union,
  21. # shapely.coverage_union is tested separately
  22. )
  23. REDUCE_SET_OPERATIONS = (
  24. (shapely.intersection_all, shapely.intersection),
  25. (shapely.symmetric_difference_all, shapely.symmetric_difference),
  26. (shapely.union_all, shapely.union),
  27. # shapely.coverage_union_all, shapely.coverage_union) is tested separately
  28. )
  29. # operations that support fixed precision
  30. REDUCE_SET_OPERATIONS_PREC = ((shapely.union_all, shapely.union),)
  31. if shapely.geos_version >= (3, 12, 0):
  32. SET_OPERATIONS += (shapely.disjoint_subset_union,)
  33. REDUCE_SET_OPERATIONS += (
  34. (shapely.disjoint_subset_union_all, shapely.disjoint_subset_union),
  35. )
  36. reduce_test_data = [
  37. shapely.box(0, 0, 5, 5),
  38. shapely.box(2, 2, 7, 7),
  39. shapely.box(4, 4, 9, 9),
  40. shapely.box(5, 5, 10, 10),
  41. ]
  42. non_polygon_types = np.array(all_types)[
  43. ~shapely.is_empty(all_types) & (shapely.get_dimensions(all_types) != 2)
  44. ]
  45. @pytest.mark.parametrize("a", all_types)
  46. @pytest.mark.parametrize("func", SET_OPERATIONS)
  47. def test_set_operation_array(a, func):
  48. if (
  49. func is shapely.difference
  50. and a.geom_type == "GeometryCollection"
  51. and shapely.get_num_geometries(a) == 2
  52. and shapely.geos_version == (3, 9, 5)
  53. ):
  54. pytest.xfail("GEOS 3.9.5 crashes with mixed collection")
  55. actual = func(a, point)
  56. assert isinstance(actual, Geometry)
  57. actual = func([a, a], point)
  58. assert actual.shape == (2,)
  59. assert isinstance(actual[0], Geometry)
  60. @pytest.mark.parametrize("func", SET_OPERATIONS)
  61. def test_set_operation_prec_nonscalar_grid_size(func):
  62. if func is shapely.disjoint_subset_union:
  63. pytest.skip("disjoint_subset_union does not support grid_size")
  64. with pytest.raises(
  65. ValueError, match="grid_size parameter only accepts scalar values"
  66. ):
  67. func(point, point, grid_size=[1])
  68. @pytest.mark.parametrize("a", all_single_types)
  69. @pytest.mark.parametrize("func", SET_OPERATIONS)
  70. @pytest.mark.parametrize("grid_size", [0, 1, 2])
  71. def test_set_operation_prec_array(a, func, grid_size):
  72. if func is shapely.disjoint_subset_union:
  73. pytest.skip("disjoint_subset_union does not support grid_size")
  74. actual = func([a, a], point, grid_size=grid_size)
  75. assert actual.shape == (2,)
  76. assert isinstance(actual[0], Geometry)
  77. # results should match the operation when the precision is previously set
  78. # to same grid_size
  79. b = shapely.set_precision(a, grid_size=grid_size)
  80. point2 = shapely.set_precision(point, grid_size=grid_size)
  81. expected = func([b, b], point2)
  82. assert shapely.equals(shapely.normalize(actual), shapely.normalize(expected)).all()
  83. @pytest.mark.parametrize("n", range(1, 5))
  84. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  85. def test_set_operation_reduce_1dim(n, func, related_func):
  86. actual = func(reduce_test_data[:n])
  87. # perform the reduction in a python loop and compare
  88. expected = reduce_test_data[0]
  89. for i in range(1, n):
  90. expected = related_func(expected, reduce_test_data[i])
  91. assert shapely.equals(actual, expected)
  92. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  93. def test_set_operation_reduce_single_geom(func, related_func):
  94. geom = shapely.Point(1, 1)
  95. actual = func([geom, None, None])
  96. assert shapely.equals(actual, geom)
  97. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  98. def test_set_operation_reduce_axis(func, related_func):
  99. data = [[point] * 2] * 3 # shape = (3, 2)
  100. actual = func(data, axis=None) # default
  101. assert isinstance(actual, Geometry) # scalar output
  102. actual = func(data, axis=0)
  103. assert actual.shape == (2,)
  104. actual = func(data, axis=1)
  105. assert actual.shape == (3,)
  106. actual = func(data, axis=-1)
  107. assert actual.shape == (3,)
  108. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  109. def test_set_operation_reduce_empty(func, related_func):
  110. assert func(np.empty((0,), dtype=object)) == empty
  111. arr_empty_2D = np.empty((0, 2), dtype=object)
  112. assert func(arr_empty_2D) == empty
  113. assert func(arr_empty_2D, axis=0).tolist() == [empty] * 2
  114. assert func(arr_empty_2D, axis=1).tolist() == []
  115. @pytest.mark.parametrize("none_position", range(3))
  116. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  117. def test_set_operation_reduce_one_none(func, related_func, none_position):
  118. # API change: before, intersection_all and symmetric_difference_all returned
  119. # None if any input geometry was None.
  120. # The new behaviour is to ignore None values.
  121. test_data = reduce_test_data[:2]
  122. test_data.insert(none_position, None)
  123. actual = func(test_data)
  124. expected = related_func(reduce_test_data[0], reduce_test_data[1])
  125. assert_geometries_equal(actual, expected)
  126. @pytest.mark.parametrize("none_position", range(3))
  127. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  128. def test_set_operation_reduce_two_none(func, related_func, none_position):
  129. test_data = reduce_test_data[:2]
  130. test_data.insert(none_position, None)
  131. test_data.insert(none_position, None)
  132. actual = func(test_data)
  133. expected = related_func(reduce_test_data[0], reduce_test_data[1])
  134. assert_geometries_equal(actual, expected)
  135. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  136. def test_set_operation_reduce_some_none_len2(func, related_func):
  137. # in a previous implementation, this would take a different code path
  138. # and return wrong result
  139. assert func([empty, None]) == empty
  140. @pytest.mark.parametrize("n", range(1, 3))
  141. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  142. def test_set_operation_reduce_all_none(n, func, related_func):
  143. assert_geometries_equal(func([None] * n), GeometryCollection([]))
  144. @pytest.mark.parametrize("n", range(1, 3))
  145. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS)
  146. def test_set_operation_reduce_all_none_arr(n, func, related_func):
  147. assert func([[None] * n] * 2, axis=1).tolist() == [empty, empty]
  148. assert func([[None] * 2] * n, axis=0).tolist() == [empty, empty]
  149. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  150. def test_set_operation_prec_reduce_nonscalar_grid_size(func, related_func):
  151. with pytest.raises(
  152. ValueError, match="grid_size parameter only accepts scalar values"
  153. ):
  154. func([point, point], grid_size=[1])
  155. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  156. def test_set_operation_prec_reduce_grid_size_nan(func, related_func):
  157. actual = func([point, point], grid_size=np.nan)
  158. assert actual is None
  159. @pytest.mark.parametrize("n", range(1, 5))
  160. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  161. @pytest.mark.parametrize("grid_size", [0, 1])
  162. def test_set_operation_prec_reduce_1dim(n, func, related_func, grid_size):
  163. actual = func(reduce_test_data[:n], grid_size=grid_size)
  164. # perform the reduction in a python loop and compare
  165. expected = reduce_test_data[0]
  166. for i in range(1, n):
  167. expected = related_func(expected, reduce_test_data[i], grid_size=grid_size)
  168. assert shapely.equals(actual, expected)
  169. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  170. def test_set_operation_prec_reduce_axis(func, related_func):
  171. data = [[point] * 2] * 3 # shape = (3, 2)
  172. actual = func(data, grid_size=1, axis=None) # default
  173. assert isinstance(actual, Geometry) # scalar output
  174. actual = func(data, grid_size=1, axis=0)
  175. assert actual.shape == (2,)
  176. actual = func(data, grid_size=1, axis=1)
  177. assert actual.shape == (3,)
  178. actual = func(data, grid_size=1, axis=-1)
  179. assert actual.shape == (3,)
  180. @pytest.mark.parametrize("none_position", range(3))
  181. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  182. def test_set_operation_prec_reduce_one_none(func, related_func, none_position):
  183. test_data = reduce_test_data[:2]
  184. test_data.insert(none_position, None)
  185. actual = func(test_data, grid_size=1)
  186. expected = related_func(reduce_test_data[0], reduce_test_data[1], grid_size=1)
  187. assert_geometries_equal(actual, expected)
  188. @pytest.mark.parametrize("none_position", range(3))
  189. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  190. def test_set_operation_prec_reduce_two_none(func, related_func, none_position):
  191. test_data = reduce_test_data[:2]
  192. test_data.insert(none_position, None)
  193. test_data.insert(none_position, None)
  194. actual = func(test_data, grid_size=1)
  195. expected = related_func(reduce_test_data[0], reduce_test_data[1], grid_size=1)
  196. assert_geometries_equal(actual, expected)
  197. @pytest.mark.parametrize("n", range(1, 3))
  198. @pytest.mark.parametrize("func, related_func", REDUCE_SET_OPERATIONS_PREC)
  199. def test_set_operation_prec_reduce_all_none(n, func, related_func):
  200. assert_geometries_equal(func([None] * n, grid_size=1), GeometryCollection([]))
  201. @pytest.mark.parametrize("n", range(1, 4))
  202. def test_coverage_union_reduce_1dim(n):
  203. """
  204. This is tested separately from other set operations as it expects only
  205. non-overlapping polygons
  206. """
  207. test_data = [
  208. shapely.box(0, 0, 1, 1),
  209. shapely.box(1, 0, 2, 1),
  210. shapely.box(2, 0, 3, 1),
  211. ]
  212. actual = shapely.coverage_union_all(test_data[:n])
  213. # perform the reduction in a python loop and compare
  214. expected = test_data[0]
  215. for i in range(1, n):
  216. expected = shapely.coverage_union(expected, test_data[i])
  217. assert_geometries_equal(actual, expected, normalize=True)
  218. def test_coverage_union_reduce_axis():
  219. # shape = (3, 2), all polygons - none of them overlapping
  220. data = [[shapely.box(i, j, i + 1, j + 1) for i in range(2)] for j in range(3)]
  221. actual = shapely.coverage_union_all(data, axis=None) # default
  222. assert isinstance(actual, Geometry)
  223. actual = shapely.coverage_union_all(data, axis=0)
  224. assert actual.shape == (2,)
  225. actual = shapely.coverage_union_all(data, axis=1)
  226. assert actual.shape == (3,)
  227. actual = shapely.coverage_union_all(data, axis=-1)
  228. assert actual.shape == (3,)
  229. def test_coverage_union_overlapping_inputs():
  230. polygon = Polygon([(1, 1), (1, 0), (0, 0), (0, 1), (1, 1)])
  231. other = Polygon([(1, 0), (0.9, 1), (2, 1), (2, 0), (1, 0)])
  232. if shapely.geos_version >= (3, 14, 0):
  233. # Overlapping polygons raise an error again
  234. with pytest.raises(shapely.GEOSException, match="TopologyException"):
  235. shapely.coverage_union(polygon, other)
  236. elif shapely.geos_version >= (3, 12, 0):
  237. # Return mostly unchanged output
  238. result = shapely.coverage_union(polygon, other)
  239. expected = shapely.multipolygons([polygon, other])
  240. assert_geometries_equal(result, expected, normalize=True)
  241. else:
  242. # Overlapping polygons raise an error
  243. with pytest.raises(
  244. shapely.GEOSException,
  245. match="CoverageUnion cannot process incorrectly noded inputs.",
  246. ):
  247. shapely.coverage_union(polygon, other)
  248. @pytest.mark.parametrize(
  249. "geom_1, geom_2",
  250. # All possible polygon, non_polygon combinations
  251. [[polygon, non_polygon] for non_polygon in non_polygon_types]
  252. # All possible non_polygon, non_polygon combinations
  253. + [
  254. [non_polygon_1, non_polygon_2]
  255. for non_polygon_1 in non_polygon_types
  256. for non_polygon_2 in non_polygon_types
  257. ],
  258. )
  259. def test_coverage_union_non_polygon_inputs(geom_1, geom_2):
  260. if shapely.geos_version >= (3, 12, 0):
  261. def effective_geom_types(geom):
  262. if hasattr(geom, "geoms") and not geom.is_empty:
  263. gts = set()
  264. for part in geom.geoms:
  265. gts |= effective_geom_types(part)
  266. return gts
  267. return {geom.geom_type.lstrip("Multi").replace("LinearRing", "LineString")}
  268. geom_types_1 = effective_geom_types(geom_1)
  269. geom_types_2 = effective_geom_types(geom_2)
  270. if len(geom_types_1) == 1 and geom_types_1 == geom_types_2:
  271. with ignore_invalid():
  272. # these show "invalid value encountered in coverage_union"
  273. result = shapely.coverage_union(geom_1, geom_2)
  274. assert geom_types_1 == effective_geom_types(result)
  275. else:
  276. with pytest.raises(
  277. shapely.GEOSException, match="Overlay input is mixed-dimension"
  278. ):
  279. shapely.coverage_union(geom_1, geom_2)
  280. else:
  281. # Non polygon geometries raise an error
  282. with pytest.raises(
  283. shapely.GEOSException, match="Unhandled geometry type in CoverageUnion."
  284. ):
  285. shapely.coverage_union(geom_1, geom_2)
  286. @pytest.mark.parametrize(
  287. "geom,grid_size,expected",
  288. [
  289. # floating point precision, expect no change
  290. (
  291. [shapely.box(0.1, 0.1, 5, 5), shapely.box(0, 0.2, 5.1, 10)],
  292. 0,
  293. Polygon(
  294. (
  295. (0, 0.2),
  296. (0, 10),
  297. (5.1, 10),
  298. (5.1, 0.2),
  299. (5, 0.2),
  300. (5, 0.1),
  301. (0.1, 0.1),
  302. (0.1, 0.2),
  303. (0, 0.2),
  304. )
  305. ),
  306. ),
  307. # grid_size is at effective precision, expect no change
  308. (
  309. [shapely.box(0.1, 0.1, 5, 5), shapely.box(0, 0.2, 5.1, 10)],
  310. 0.1,
  311. Polygon(
  312. (
  313. (0, 0.2),
  314. (0, 10),
  315. (5.1, 10),
  316. (5.1, 0.2),
  317. (5, 0.2),
  318. (5, 0.1),
  319. (0.1, 0.1),
  320. (0.1, 0.2),
  321. (0, 0.2),
  322. )
  323. ),
  324. ),
  325. # grid_size forces rounding to nearest integer
  326. (
  327. [shapely.box(0.1, 0.1, 5, 5), shapely.box(0, 0.2, 5.1, 10)],
  328. 1,
  329. Polygon([(0, 5), (0, 10), (5, 10), (5, 5), (5, 0), (0, 0), (0, 5)]),
  330. ),
  331. # grid_size much larger than effective precision causes rounding to nearest
  332. # multiple of 10
  333. (
  334. [shapely.box(0.1, 0.1, 5, 5), shapely.box(0, 0.2, 5.1, 10)],
  335. 10,
  336. Polygon([(0, 10), (10, 10), (10, 0), (0, 0), (0, 10)]),
  337. ),
  338. # grid_size is so large that polygons collapse to empty
  339. (
  340. [shapely.box(0.1, 0.1, 5, 5), shapely.box(0, 0.2, 5.1, 10)],
  341. 100,
  342. Polygon(),
  343. ),
  344. ],
  345. )
  346. def test_union_all_prec(geom, grid_size, expected):
  347. actual = shapely.union_all(geom, grid_size=grid_size)
  348. assert shapely.equals(actual, expected)
  349. def test_uary_union_alias():
  350. geoms = [shapely.box(0.1, 0.1, 5, 5), shapely.box(0, 0.2, 5.1, 10)]
  351. actual = shapely.unary_union(geoms, grid_size=1)
  352. expected = shapely.union_all(geoms, grid_size=1)
  353. assert shapely.equals(actual, expected)
  354. def test_difference_deprecate_positional():
  355. with pytest.deprecated_call(
  356. match="positional argument `grid_size` for `difference` is deprecated"
  357. ):
  358. shapely.difference(point, point, None)
  359. def test_intersection_deprecate_positional():
  360. with pytest.deprecated_call(
  361. match="positional argument `grid_size` for `intersection` is deprecated"
  362. ):
  363. shapely.intersection(point, point, None)
  364. def test_intersection_all_deprecate_positional():
  365. with pytest.deprecated_call(
  366. match="positional argument `axis` for `intersection_all` is deprecated"
  367. ):
  368. shapely.intersection_all([point, point], None)
  369. def test_symmetric_difference_deprecate_positional():
  370. with pytest.deprecated_call(
  371. match="positional argument `grid_size` for `symmetric_difference` is deprecated"
  372. ):
  373. shapely.symmetric_difference(point, point, None)
  374. def test_symmetric_difference_all_deprecate_positional():
  375. with pytest.deprecated_call(
  376. match="positional argument `axis` for `symmetric_difference_all` is deprecated"
  377. ):
  378. shapely.symmetric_difference_all([point, point], None)
  379. def test_union_deprecate_positional():
  380. with pytest.deprecated_call(
  381. match="positional argument `grid_size` for `union` is deprecated"
  382. ):
  383. shapely.union(point, point, None)
  384. def test_union_all_deprecate_positional():
  385. with pytest.deprecated_call(
  386. match="positional argument `grid_size` for `union_all` is deprecated"
  387. ):
  388. shapely.union_all([point, point], None)
  389. with pytest.deprecated_call(
  390. match="positional arguments `grid_size` and `axis` for `union_all` "
  391. "are deprecated"
  392. ):
  393. shapely.union_all([point, point], None, None)
  394. def test_coverage_union_all_deprecate_positional():
  395. data = [shapely.box(0, 0, 1, 1), shapely.box(1, 0, 2, 1)]
  396. with pytest.deprecated_call(
  397. match="positional argument `axis` for `coverage_union_all` is deprecated"
  398. ):
  399. shapely.coverage_union_all(data, None)