test_strtree.py 72 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964
  1. import itertools
  2. import math
  3. import pickle
  4. import subprocess
  5. import sys
  6. from concurrent.futures import ThreadPoolExecutor
  7. import numpy as np
  8. import pytest
  9. from numpy.testing import assert_array_equal
  10. import shapely
  11. from shapely import LineString, MultiPoint, Point, STRtree, box, geos_version
  12. from shapely.errors import UnsupportedGEOSVersionError
  13. from shapely.testing import assert_geometries_equal
  14. from shapely.tests.common import (
  15. empty,
  16. empty_line_string,
  17. empty_point,
  18. ignore_invalid,
  19. point,
  20. )
  21. # the distance between 2 points spaced at whole numbers along a diagonal
  22. HALF_UNIT_DIAG = math.sqrt(2) / 2
  23. EPS = 1e-9
  24. @pytest.fixture(scope="session")
  25. def tree():
  26. geoms = shapely.points(np.arange(10), np.arange(10))
  27. yield STRtree(geoms)
  28. @pytest.fixture(scope="session")
  29. def line_tree():
  30. x = np.arange(10)
  31. y = np.arange(10)
  32. offset = 1
  33. geoms = shapely.linestrings(np.array([[x, x + offset], [y, y + offset]]).T)
  34. yield STRtree(geoms)
  35. @pytest.fixture(scope="session")
  36. def poly_tree():
  37. # create buffers so that midpoint between two buffers intersects
  38. # each buffer. NOTE: add EPS to help mitigate rounding errors at midpoint.
  39. geoms = shapely.buffer(
  40. shapely.points(np.arange(10), np.arange(10)), HALF_UNIT_DIAG + EPS, quad_segs=32
  41. )
  42. yield STRtree(geoms)
  43. @pytest.mark.parametrize(
  44. "geometry,count, hits",
  45. [
  46. # Empty array produces empty tree
  47. ([], 0, 0),
  48. ([point], 1, 1),
  49. # None geometries are ignored when creating tree
  50. ([None], 0, 0),
  51. ([point, None], 1, 1),
  52. # empty geometries are ignored when creating tree
  53. ([empty, empty_point, empty_line_string], 0, 0),
  54. # only the valid geometry should have a hit
  55. ([empty, point, empty_point, empty_line_string], 1, 1),
  56. ],
  57. )
  58. def test_init(geometry, count, hits):
  59. tree = STRtree(geometry)
  60. assert len(tree) == count
  61. assert tree.query(box(0, 0, 100, 100)).size == hits
  62. def test_init_with_invalid_geometry():
  63. with pytest.raises(TypeError):
  64. STRtree(["Not a geometry"])
  65. def test_references():
  66. point1 = Point()
  67. point2 = Point(0, 1)
  68. geoms = [point1, point2]
  69. tree = STRtree(geoms)
  70. point1 = None
  71. point2 = None
  72. import gc
  73. gc.collect()
  74. # query after freeing geometries does not lead to segfault
  75. assert tree.query(box(0, 0, 1, 1)).tolist() == [1]
  76. def test_flush_geometries():
  77. arr = shapely.points(np.arange(10), np.arange(10))
  78. tree = STRtree(arr)
  79. # Dereference geometries
  80. arr[:] = None
  81. import gc
  82. gc.collect()
  83. # Still it does not lead to a segfault
  84. tree.query(point)
  85. def test_geometries_property():
  86. arr = np.array([point])
  87. tree = STRtree(arr)
  88. assert_geometries_equal(arr, tree.geometries)
  89. # modifying elements of input should not modify tree.geometries
  90. arr[0] = shapely.Point(0, 0)
  91. assert_geometries_equal(point, tree.geometries[0])
  92. def test_pickle_persistence(tmp_path):
  93. # write the pickeled tree to another process; the process should not crash
  94. tree = STRtree([Point(i, i).buffer(0.1) for i in range(3)])
  95. pickled_strtree = pickle.dumps(tree)
  96. unpickle_script = """
  97. import pickle
  98. import sys
  99. from shapely import Point
  100. pickled_strtree = sys.stdin.buffer.read()
  101. print("received pickled strtree:", repr(pickled_strtree))
  102. tree = pickle.loads(pickled_strtree)
  103. tree.query(Point(0, 0))
  104. tree.nearest(Point(0, 0))
  105. print("done")
  106. """
  107. filename = tmp_path / "unpickle-strtree.py"
  108. with open(filename, "w") as out:
  109. out.write(unpickle_script)
  110. proc = subprocess.Popen(
  111. [sys.executable, str(filename)],
  112. stdin=subprocess.PIPE,
  113. )
  114. proc.communicate(input=pickled_strtree)
  115. proc.wait()
  116. assert proc.returncode == 0
  117. @pytest.mark.parametrize(
  118. "geometry",
  119. [
  120. "I am not a geometry",
  121. ["I am not a geometry"],
  122. [Point(0, 0), "still not a geometry"],
  123. [[], "in a mixed array", 1],
  124. ],
  125. )
  126. @pytest.mark.filterwarnings("ignore:Creating an ndarray from ragged nested sequences:")
  127. def test_query_invalid_geometry(tree, geometry):
  128. with pytest.raises((TypeError, ValueError)):
  129. tree.query(geometry)
  130. def test_query_invalid_dimension(tree):
  131. with pytest.raises(TypeError, match="Array should be one dimensional"):
  132. tree.query([[Point(0.5, 0.5)]])
  133. @pytest.mark.parametrize(
  134. "tree_geometry, geometry,expected",
  135. [
  136. # Empty tree returns no results
  137. ([], point, []),
  138. ([], [point], [[], []]),
  139. ([], None, []),
  140. ([], [None], [[], []]),
  141. # Tree with only None returns no results
  142. ([None], point, []),
  143. ([None], [point], [[], []]),
  144. ([None], None, []),
  145. ([None], [None], [[], []]),
  146. # querying with None returns no results
  147. ([point], None, []),
  148. ([point], [None], [[], []]),
  149. # Empty is included in the tree, but ignored when querying the tree
  150. ([empty], empty, []),
  151. ([empty], [empty], [[], []]),
  152. ([empty], point, []),
  153. ([empty], [point], [[], []]),
  154. ([point, empty], empty, []),
  155. ([point, empty], [empty], [[], []]),
  156. # None and empty are ignored in the tree, but the index of the valid
  157. # geometry should be retained.
  158. ([None, point], box(0, 0, 10, 10), [1]),
  159. ([None, point], [box(0, 0, 10, 10)], [[0], [1]]),
  160. ([None, empty, point], box(0, 0, 10, 10), [2]),
  161. ([point, None, point], box(0, 0, 10, 10), [0, 2]),
  162. ([point, None, point], [box(0, 0, 10, 10)], [[0, 0], [0, 2]]),
  163. # Only the non-empty query geometry gets hits
  164. ([empty, point], [empty, point], [[1], [1]]),
  165. (
  166. [empty, empty_point, empty_line_string, point],
  167. [empty, empty_point, empty_line_string, point],
  168. [[3], [3]],
  169. ),
  170. ],
  171. )
  172. def test_query_with_none_and_empty(tree_geometry, geometry, expected):
  173. tree = STRtree(tree_geometry)
  174. assert_array_equal(tree.query(geometry), expected)
  175. @pytest.mark.parametrize(
  176. "geometry,expected",
  177. [
  178. # points do not intersect
  179. (Point(0.5, 0.5), []),
  180. ([Point(0.5, 0.5)], [[], []]),
  181. # points intersect
  182. (Point(1, 1), [1]),
  183. ([Point(1, 1)], [[0], [1]]),
  184. # first and last points intersect
  185. (
  186. [Point(1, 1), Point(-1, -1), Point(2, 2)],
  187. [[0, 2], [1, 2]],
  188. ),
  189. # box contains points
  190. (box(0, 0, 1, 1), [0, 1]),
  191. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  192. # bigger box contains more points
  193. (box(5, 5, 15, 15), [5, 6, 7, 8, 9]),
  194. ([box(5, 5, 15, 15)], [[0, 0, 0, 0, 0], [5, 6, 7, 8, 9]]),
  195. # first and last boxes contains points
  196. (
  197. [box(0, 0, 1, 1), box(100, 100, 110, 110), box(5, 5, 15, 15)],
  198. [[0, 0, 2, 2, 2, 2, 2], [0, 1, 5, 6, 7, 8, 9]],
  199. ),
  200. # envelope of buffer contains points
  201. (shapely.buffer(Point(3, 3), 1), [2, 3, 4]),
  202. ([shapely.buffer(Point(3, 3), 1)], [[0, 0, 0], [2, 3, 4]]),
  203. # envelope of points contains points
  204. (MultiPoint([[5, 7], [7, 5]]), [5, 6, 7]),
  205. ([MultiPoint([[5, 7], [7, 5]])], [[0, 0, 0], [5, 6, 7]]),
  206. ],
  207. )
  208. def test_query_points(tree, geometry, expected):
  209. assert_array_equal(tree.query(geometry), expected)
  210. @pytest.mark.parametrize(
  211. "geometry,expected",
  212. [
  213. # point intersects first line
  214. (Point(0, 0), [0]),
  215. ([Point(0, 0)], [[0], [0]]),
  216. (Point(0.5, 0.5), [0]),
  217. ([Point(0.5, 0.5)], [[0], [0]]),
  218. # point within envelope of first line
  219. (Point(0, 0.5), [0]),
  220. ([Point(0, 0.5)], [[0], [0]]),
  221. # point at shared vertex between 2 lines
  222. (Point(1, 1), [0, 1]),
  223. ([Point(1, 1)], [[0, 0], [0, 1]]),
  224. # box overlaps envelope of first 2 lines (touches edge of 1)
  225. (box(0, 0, 1, 1), [0, 1]),
  226. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  227. # envelope of buffer overlaps envelope of 2 lines
  228. (shapely.buffer(Point(3, 3), 0.5), [2, 3]),
  229. ([shapely.buffer(Point(3, 3), 0.5)], [[0, 0], [2, 3]]),
  230. # envelope of points overlaps 5 lines (touches edge of 2 envelopes)
  231. (MultiPoint([[5, 7], [7, 5]]), [4, 5, 6, 7]),
  232. ([MultiPoint([[5, 7], [7, 5]])], [[0, 0, 0, 0], [4, 5, 6, 7]]),
  233. ],
  234. )
  235. def test_query_lines(line_tree, geometry, expected):
  236. assert_array_equal(line_tree.query(geometry), expected)
  237. @pytest.mark.parametrize(
  238. "geometry,expected",
  239. [
  240. # point intersects edge of envelopes of 2 polygons
  241. (Point(0.5, 0.5), [0, 1]),
  242. ([Point(0.5, 0.5)], [[0, 0], [0, 1]]),
  243. # point intersects single polygon
  244. (Point(1, 1), [1]),
  245. ([Point(1, 1)], [[0], [1]]),
  246. # box overlaps envelope of 2 polygons
  247. (box(0, 0, 1, 1), [0, 1]),
  248. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  249. # larger box overlaps envelope of 3 polygons
  250. (box(0, 0, 1.5, 1.5), [0, 1, 2]),
  251. ([box(0, 0, 1.5, 1.5)], [[0, 0, 0], [0, 1, 2]]),
  252. # first and last boxes overlap envelope of 2 polyons
  253. (
  254. [box(0, 0, 1, 1), box(100, 100, 110, 110), box(2, 2, 3, 3)],
  255. [[0, 0, 2, 2], [0, 1, 2, 3]],
  256. ),
  257. # envelope of buffer overlaps envelope of 3 polygons
  258. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [2, 3, 4]),
  259. (
  260. [shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)],
  261. [[0, 0, 0], [2, 3, 4]],
  262. ),
  263. # envelope of larger buffer overlaps envelope of 6 polygons
  264. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [1, 2, 3, 4, 5]),
  265. (
  266. [shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)],
  267. [[0, 0, 0, 0, 0], [1, 2, 3, 4, 5]],
  268. ),
  269. # envelope of points overlaps 3 polygons
  270. (MultiPoint([[5, 7], [7, 5]]), [5, 6, 7]),
  271. ([MultiPoint([[5, 7], [7, 5]])], [[0, 0, 0], [5, 6, 7]]),
  272. ],
  273. )
  274. def test_query_polygons(poly_tree, geometry, expected):
  275. assert_array_equal(poly_tree.query(geometry), expected)
  276. @pytest.mark.parametrize(
  277. "predicate",
  278. [
  279. "bad_predicate",
  280. # disjoint is a valid GEOS binary predicate, but not supported for query
  281. "disjoint",
  282. ],
  283. )
  284. def test_query_invalid_predicate(tree, predicate):
  285. with pytest.raises(ValueError, match="is not a valid option"):
  286. tree.query(Point(1, 1), predicate=predicate)
  287. @pytest.mark.parametrize(
  288. "predicate,expected",
  289. [
  290. ("intersects", [0, 1, 2]),
  291. ("within", []),
  292. ("contains", [1]),
  293. ("overlaps", []),
  294. ("crosses", []),
  295. ("covers", [0, 1, 2]),
  296. ("covered_by", []),
  297. ("contains_properly", [1]),
  298. ],
  299. )
  300. def test_query_prepared_inputs(tree, predicate, expected):
  301. geom = box(0, 0, 2, 2)
  302. shapely.prepare(geom)
  303. assert_array_equal(tree.query(geom, predicate=predicate), expected)
  304. def test_query_with_partially_prepared_inputs(tree):
  305. geom = np.array([box(0, 0, 1, 1), box(3, 3, 5, 5)])
  306. expected = tree.query(geom, predicate="intersects")
  307. # test with array of partially prepared geometries
  308. shapely.prepare(geom[0])
  309. assert_array_equal(expected, tree.query(geom, predicate="intersects"))
  310. @pytest.mark.parametrize(
  311. "predicate,expected",
  312. [
  313. pytest.param(
  314. "intersects",
  315. [1],
  316. marks=pytest.mark.xfail(geos_version < (3, 13, 0), reason="GEOS < 3.13"),
  317. ),
  318. ("within", []),
  319. ("contains", []),
  320. ("overlaps", []),
  321. ("crosses", [1]),
  322. ("touches", []),
  323. ("covers", []),
  324. ("covered_by", []),
  325. ("contains_properly", []),
  326. ],
  327. )
  328. def test_query_predicate_errors(tree, predicate, expected):
  329. with ignore_invalid():
  330. line_nan = shapely.linestrings([1, 1], [1, float("nan")])
  331. if geos_version < (3, 13, 0):
  332. with pytest.raises(shapely.GEOSException):
  333. tree.query(line_nan, predicate=predicate)
  334. else:
  335. assert_array_equal(tree.query(line_nan, predicate=predicate), expected)
  336. ### predicate == 'intersects'
  337. @pytest.mark.parametrize(
  338. "geometry,expected",
  339. [
  340. # points do not intersect
  341. (Point(0.5, 0.5), []),
  342. ([Point(0.5, 0.5)], [[], []]),
  343. # points intersect
  344. (Point(1, 1), [1]),
  345. ([Point(1, 1)], [[0], [1]]),
  346. # box contains points
  347. (box(3, 3, 6, 6), [3, 4, 5, 6]),
  348. ([box(3, 3, 6, 6)], [[0, 0, 0, 0], [3, 4, 5, 6]]),
  349. # first and last boxes contain points
  350. (
  351. [box(0, 0, 1, 1), box(100, 100, 110, 110), box(3, 3, 6, 6)],
  352. [[0, 0, 2, 2, 2, 2], [0, 1, 3, 4, 5, 6]],
  353. ),
  354. # envelope of buffer contains more points than intersect buffer
  355. # due to diagonal distance
  356. (shapely.buffer(Point(3, 3), 1), [3]),
  357. ([shapely.buffer(Point(3, 3), 1)], [[0], [3]]),
  358. # envelope of buffer with 1/2 distance between points should intersect
  359. # same points as envelope
  360. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [2, 3, 4]),
  361. (
  362. [shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)],
  363. [[0, 0, 0], [2, 3, 4]],
  364. ),
  365. # multipoints intersect
  366. (
  367. MultiPoint([[5, 5], [7, 7]]),
  368. [5, 7],
  369. ),
  370. (
  371. [MultiPoint([[5, 5], [7, 7]])],
  372. [[0, 0], [5, 7]],
  373. ),
  374. # envelope of points contains points, but points do not intersect
  375. (MultiPoint([[5, 7], [7, 5]]), []),
  376. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  377. # only one point of multipoint intersects
  378. (
  379. MultiPoint([[5, 7], [7, 7]]),
  380. [7],
  381. ),
  382. (
  383. [MultiPoint([[5, 7], [7, 7]])],
  384. [[0], [7]],
  385. ),
  386. ],
  387. )
  388. def test_query_intersects_points(tree, geometry, expected):
  389. assert_array_equal(tree.query(geometry, predicate="intersects"), expected)
  390. @pytest.mark.parametrize(
  391. "geometry,expected",
  392. [
  393. # point intersects first line
  394. (Point(0, 0), [0]),
  395. ([Point(0, 0)], [[0], [0]]),
  396. (Point(0.5, 0.5), [0]),
  397. ([Point(0.5, 0.5)], [[0], [0]]),
  398. # point within envelope of first line but does not intersect
  399. (Point(0, 0.5), []),
  400. ([Point(0, 0.5)], [[], []]),
  401. # point at shared vertex between 2 lines
  402. (Point(1, 1), [0, 1]),
  403. ([Point(1, 1)], [[0, 0], [0, 1]]),
  404. # box overlaps envelope of first 2 lines (touches edge of 1)
  405. (box(0, 0, 1, 1), [0, 1]),
  406. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  407. # first and last boxes overlap multiple lines each
  408. (
  409. [box(0, 0, 1, 1), box(100, 100, 110, 110), box(2, 2, 3, 3)],
  410. [[0, 0, 2, 2, 2], [0, 1, 1, 2, 3]],
  411. ),
  412. # buffer intersects 2 lines
  413. (shapely.buffer(Point(3, 3), 0.5), [2, 3]),
  414. ([shapely.buffer(Point(3, 3), 0.5)], [[0, 0], [2, 3]]),
  415. # buffer intersects midpoint of line at tangent
  416. (shapely.buffer(Point(2, 1), HALF_UNIT_DIAG), [1]),
  417. ([shapely.buffer(Point(2, 1), HALF_UNIT_DIAG)], [[0], [1]]),
  418. # envelope of points overlaps lines but intersects none
  419. (MultiPoint([[5, 7], [7, 5]]), []),
  420. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  421. # only one point of multipoint intersects
  422. (MultiPoint([[5, 7], [7, 7]]), [6, 7]),
  423. ([MultiPoint([[5, 7], [7, 7]])], [[0, 0], [6, 7]]),
  424. ],
  425. )
  426. def test_query_intersects_lines(line_tree, geometry, expected):
  427. assert_array_equal(line_tree.query(geometry, predicate="intersects"), expected)
  428. @pytest.mark.parametrize(
  429. "geometry,expected",
  430. [
  431. # point within first polygon
  432. (Point(0, 0.5), [0]),
  433. ([Point(0, 0.5)], [[0], [0]]),
  434. (Point(0.5, 0), [0]),
  435. ([Point(0.5, 0)], [[0], [0]]),
  436. # midpoint between two polygons intersects both
  437. (Point(0.5, 0.5), [0, 1]),
  438. ([Point(0.5, 0.5)], [[0, 0], [0, 1]]),
  439. # point intersects single polygon
  440. (Point(1, 1), [1]),
  441. ([Point(1, 1)], [[0], [1]]),
  442. # box overlaps envelope of 2 polygons
  443. (box(0, 0, 1, 1), [0, 1]),
  444. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  445. # larger box intersects 3 polygons
  446. (box(0, 0, 1.5, 1.5), [0, 1, 2]),
  447. ([box(0, 0, 1.5, 1.5)], [[0, 0, 0], [0, 1, 2]]),
  448. # first and last boxes overlap
  449. (
  450. [box(0, 0, 1, 1), box(100, 100, 110, 110), box(2, 2, 3, 3)],
  451. [[0, 0, 2, 2], [0, 1, 2, 3]],
  452. ),
  453. # buffer overlaps 3 polygons
  454. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [2, 3, 4]),
  455. (
  456. [shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)],
  457. [[0, 0, 0], [2, 3, 4]],
  458. ),
  459. # larger buffer overlaps 6 polygons (touches midpoints)
  460. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [1, 2, 3, 4, 5]),
  461. (
  462. [shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)],
  463. [[0, 0, 0, 0, 0], [1, 2, 3, 4, 5]],
  464. ),
  465. # envelope of points overlaps polygons, but points do not intersect
  466. (MultiPoint([[5, 7], [7, 5]]), []),
  467. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  468. # only one point of multipoint within polygon
  469. (MultiPoint([[5, 7], [7, 7]]), [7]),
  470. ([MultiPoint([[5, 7], [7, 7]])], [[0], [7]]),
  471. ],
  472. )
  473. def test_query_intersects_polygons(poly_tree, geometry, expected):
  474. assert_array_equal(poly_tree.query(geometry, predicate="intersects"), expected)
  475. ### predicate == 'within'
  476. @pytest.mark.parametrize(
  477. "geometry,expected",
  478. [
  479. # points do not intersect
  480. (Point(0.5, 0.5), []),
  481. ([Point(0.5, 0.5)], [[], []]),
  482. # points intersect
  483. (Point(1, 1), [1]),
  484. ([Point(1, 1)], [[0], [1]]),
  485. # box not within points
  486. (box(3, 3, 6, 6), []),
  487. ([box(3, 3, 6, 6)], [[], []]),
  488. # envelope of buffer not within points
  489. (shapely.buffer(Point(3, 3), 1), []),
  490. ([shapely.buffer(Point(3, 3), 1)], [[], []]),
  491. # multipoints intersect but are not within points in tree
  492. (MultiPoint([[5, 5], [7, 7]]), []),
  493. ([MultiPoint([[5, 5], [7, 7]])], [[], []]),
  494. # only one point of multipoint intersects, but multipoints are not
  495. # within any points in tree
  496. (MultiPoint([[5, 7], [7, 7]]), []),
  497. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  498. # envelope of points contains points, but points do not intersect
  499. (MultiPoint([[5, 7], [7, 5]]), []),
  500. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  501. ],
  502. )
  503. def test_query_within_points(tree, geometry, expected):
  504. assert_array_equal(tree.query(geometry, predicate="within"), expected)
  505. @pytest.mark.parametrize(
  506. "geometry,expected",
  507. [
  508. # endpoint not within first line
  509. (Point(0, 0), []),
  510. ([Point(0, 0)], [[], []]),
  511. # point within first line
  512. (Point(0.5, 0.5), [0]),
  513. ([Point(0.5, 0.5)], [[0], [0]]),
  514. # point within envelope of first line but does not intersect
  515. (Point(0, 0.5), []),
  516. ([Point(0, 0.5)], [[], []]),
  517. # point at shared vertex between 2 lines (but within neither)
  518. (Point(1, 1), []),
  519. ([Point(1, 1)], [[], []]),
  520. # box not within line
  521. (box(0, 0, 1, 1), []),
  522. ([box(0, 0, 1, 1)], [[], []]),
  523. # buffer intersects 2 lines but not within either
  524. (shapely.buffer(Point(3, 3), 0.5), []),
  525. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  526. # envelope of points overlaps lines but intersects none
  527. (MultiPoint([[5, 7], [7, 5]]), []),
  528. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  529. # only one point of multipoint intersects, but both are not within line
  530. (MultiPoint([[5, 7], [7, 7]]), []),
  531. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  532. (MultiPoint([[6.5, 6.5], [7, 7]]), [6]),
  533. ([MultiPoint([[6.5, 6.5], [7, 7]])], [[0], [6]]),
  534. ],
  535. )
  536. def test_query_within_lines(line_tree, geometry, expected):
  537. assert_array_equal(line_tree.query(geometry, predicate="within"), expected)
  538. @pytest.mark.parametrize(
  539. "geometry,expected",
  540. [
  541. # point within first polygon
  542. (Point(0, 0.5), [0]),
  543. ([Point(0, 0.5)], [[0], [0]]),
  544. (Point(0.5, 0), [0]),
  545. ([Point(0.5, 0)], [[0], [0]]),
  546. # midpoint between two polygons intersects both
  547. (Point(0.5, 0.5), [0, 1]),
  548. ([Point(0.5, 0.5)], [[0, 0], [0, 1]]),
  549. # point intersects single polygon
  550. (Point(1, 1), [1]),
  551. ([Point(1, 1)], [[0], [1]]),
  552. # box overlaps envelope of 2 polygons but within neither
  553. (box(0, 0, 1, 1), []),
  554. ([box(0, 0, 1, 1)], [[], []]),
  555. # box within polygon
  556. (box(0, 0, 0.5, 0.5), [0]),
  557. ([box(0, 0, 0.5, 0.5)], [[0], [0]]),
  558. # larger box intersects 3 polygons but within none
  559. (box(0, 0, 1.5, 1.5), []),
  560. ([box(0, 0, 1.5, 1.5)], [[], []]),
  561. # buffer intersects 3 polygons but only within one
  562. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [3]),
  563. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[0], [3]]),
  564. # larger buffer overlaps 6 polygons (touches midpoints) but within none
  565. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), []),
  566. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[], []]),
  567. # envelope of points overlaps polygons, but points do not intersect
  568. (MultiPoint([[5, 7], [7, 5]]), []),
  569. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  570. # only one point of multipoint within polygon
  571. (MultiPoint([[5, 7], [7, 7]]), []),
  572. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  573. # both points in multipoint within polygon
  574. (MultiPoint([[5.25, 5.5], [5.25, 5.0]]), [5]),
  575. ([MultiPoint([[5.25, 5.5], [5.25, 5.0]])], [[0], [5]]),
  576. ],
  577. )
  578. def test_query_within_polygons(poly_tree, geometry, expected):
  579. assert_array_equal(poly_tree.query(geometry, predicate="within"), expected)
  580. ### predicate == 'contains'
  581. @pytest.mark.parametrize(
  582. "geometry,expected",
  583. [
  584. # points do not intersect
  585. (Point(0.5, 0.5), []),
  586. ([Point(0.5, 0.5)], [[], []]),
  587. # points intersect
  588. (Point(1, 1), [1]),
  589. ([Point(1, 1)], [[0], [1]]),
  590. # box contains points (2 are at edges and not contained)
  591. (box(3, 3, 6, 6), [4, 5]),
  592. ([box(3, 3, 6, 6)], [[0, 0], [4, 5]]),
  593. # envelope of buffer contains more points than within buffer
  594. # due to diagonal distance
  595. (shapely.buffer(Point(3, 3), 1), [3]),
  596. ([shapely.buffer(Point(3, 3), 1)], [[0], [3]]),
  597. # envelope of buffer with 1/2 distance between points should intersect
  598. # same points as envelope
  599. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [2, 3, 4]),
  600. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[0, 0, 0], [2, 3, 4]]),
  601. # multipoints intersect
  602. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  603. ([MultiPoint([[5, 5], [7, 7]])], [[0, 0], [5, 7]]),
  604. # envelope of points contains points, but points do not intersect
  605. (MultiPoint([[5, 7], [7, 5]]), []),
  606. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  607. # only one point of multipoint intersects
  608. (MultiPoint([[5, 7], [7, 7]]), [7]),
  609. ([MultiPoint([[5, 7], [7, 7]])], [[0], [7]]),
  610. ],
  611. )
  612. def test_query_contains_points(tree, geometry, expected):
  613. assert_array_equal(tree.query(geometry, predicate="contains"), expected)
  614. @pytest.mark.parametrize(
  615. "geometry,expected",
  616. [
  617. # point does not contain any lines (not valid relation)
  618. (Point(0, 0), []),
  619. ([Point(0, 0)], [[], []]),
  620. # box contains first line (touches edge of 1 but does not contain it)
  621. (box(0, 0, 1, 1), [0]),
  622. ([box(0, 0, 1, 1)], [[0], [0]]),
  623. # buffer intersects 2 lines but contains neither
  624. (shapely.buffer(Point(3, 3), 0.5), []),
  625. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  626. # envelope of points overlaps lines but intersects none
  627. (MultiPoint([[5, 7], [7, 5]]), []),
  628. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  629. # only one point of multipoint intersects
  630. (MultiPoint([[5, 7], [7, 7]]), []),
  631. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  632. # both points intersect but do not contain any lines (not valid relation)
  633. (MultiPoint([[5, 5], [6, 6]]), []),
  634. ([MultiPoint([[5, 5], [6, 6]])], [[], []]),
  635. ],
  636. )
  637. def test_query_contains_lines(line_tree, geometry, expected):
  638. assert_array_equal(line_tree.query(geometry, predicate="contains"), expected)
  639. @pytest.mark.parametrize(
  640. "geometry,expected",
  641. [
  642. # point does not contain any polygons (not valid relation)
  643. (Point(0, 0), []),
  644. ([Point(0, 0)], [[], []]),
  645. # box overlaps envelope of 2 polygons but contains neither
  646. (box(0, 0, 1, 1), []),
  647. ([box(0, 0, 1, 1)], [[], []]),
  648. # larger box intersects 3 polygons but contains only one
  649. (box(0, 0, 2, 2), [1]),
  650. ([box(0, 0, 2, 2)], [[0], [1]]),
  651. # buffer overlaps 3 polygons but contains none
  652. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), []),
  653. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[], []]),
  654. # larger buffer overlaps 6 polygons (touches midpoints) but contains one
  655. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [3]),
  656. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[0], [3]]),
  657. # envelope of points overlaps polygons, but points do not intersect
  658. # (not valid relation)
  659. (MultiPoint([[5, 7], [7, 5]]), []),
  660. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  661. ],
  662. )
  663. def test_query_contains_polygons(poly_tree, geometry, expected):
  664. assert_array_equal(poly_tree.query(geometry, predicate="contains"), expected)
  665. ### predicate == 'overlaps'
  666. # Overlaps only returns results where geometries are of same dimensions
  667. # and do not completely contain each other.
  668. # See: https://postgis.net/docs/ST_Overlaps.html
  669. @pytest.mark.parametrize(
  670. "geometry,expected",
  671. [
  672. # points do not intersect
  673. (Point(0.5, 0.5), []),
  674. ([Point(0.5, 0.5)], [[], []]),
  675. # points intersect but do not overlap
  676. (Point(1, 1), []),
  677. ([Point(1, 1)], [[], []]),
  678. # box overlaps points including those at edge but does not overlap
  679. # (completely contains all points)
  680. (box(3, 3, 6, 6), []),
  681. ([box(3, 3, 6, 6)], [[], []]),
  682. # envelope of buffer contains points, but does not overlap
  683. (shapely.buffer(Point(3, 3), 1), []),
  684. ([shapely.buffer(Point(3, 3), 1)], [[], []]),
  685. # multipoints intersect but do not overlap (both completely contain each other)
  686. (MultiPoint([[5, 5], [7, 7]]), []),
  687. ([MultiPoint([[5, 5], [7, 7]])], [[], []]),
  688. # envelope of points contains points in tree, but points do not intersect
  689. (MultiPoint([[5, 7], [7, 5]]), []),
  690. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  691. # only one point of multipoint intersects but does not overlap
  692. # the intersecting point from multipoint completely contains point in tree
  693. (MultiPoint([[5, 7], [7, 7]]), []),
  694. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  695. ],
  696. )
  697. def test_query_overlaps_points(tree, geometry, expected):
  698. assert_array_equal(tree.query(geometry, predicate="overlaps"), expected)
  699. @pytest.mark.parametrize(
  700. "geometry,expected",
  701. [
  702. # point intersects line but is completely contained by it
  703. (Point(0, 0), []),
  704. ([Point(0, 0)], [[], []]),
  705. # box overlaps second line (contains first line)
  706. # but of different dimensions so does not overlap
  707. (box(0, 0, 1.5, 1.5), []),
  708. ([box(0, 0, 1.5, 1.5)], [[], []]),
  709. # buffer intersects 2 lines but of different dimensions so does not overlap
  710. (shapely.buffer(Point(3, 3), 0.5), []),
  711. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  712. # envelope of points overlaps lines but intersects none
  713. (MultiPoint([[5, 7], [7, 5]]), []),
  714. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  715. # only one point of multipoint intersects
  716. (MultiPoint([[5, 7], [7, 7]]), []),
  717. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  718. # both points intersect but different dimensions
  719. (MultiPoint([[5, 5], [6, 6]]), []),
  720. ([MultiPoint([[5, 5], [6, 6]])], [[], []]),
  721. ],
  722. )
  723. def test_query_overlaps_lines(line_tree, geometry, expected):
  724. assert_array_equal(line_tree.query(geometry, predicate="overlaps"), expected)
  725. @pytest.mark.parametrize(
  726. "geometry,expected",
  727. [
  728. # point does not overlap any polygons (different dimensions)
  729. (Point(0, 0), []),
  730. ([Point(0, 0)], [[], []]),
  731. # box overlaps 2 polygons
  732. (box(0, 0, 1, 1), [0, 1]),
  733. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  734. # larger box intersects 3 polygons and contains one
  735. (box(0, 0, 2, 2), [0, 2]),
  736. ([box(0, 0, 2, 2)], [[0, 0], [0, 2]]),
  737. # buffer overlaps 3 polygons and contains 1
  738. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [2, 4]),
  739. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[0, 0], [2, 4]]),
  740. # larger buffer overlaps 6 polygons (touches midpoints) but contains one
  741. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [1, 2, 4, 5]),
  742. (
  743. [shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)],
  744. [[0, 0, 0, 0], [1, 2, 4, 5]],
  745. ),
  746. # one of two points intersects but different dimensions
  747. (MultiPoint([[5, 7], [7, 7]]), []),
  748. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  749. ],
  750. )
  751. def test_query_overlaps_polygons(poly_tree, geometry, expected):
  752. assert_array_equal(poly_tree.query(geometry, predicate="overlaps"), expected)
  753. ### predicate == 'crosses'
  754. # Only valid for certain geometry combinations
  755. # See: https://postgis.net/docs/ST_Crosses.html
  756. @pytest.mark.parametrize(
  757. "geometry,expected",
  758. [
  759. # points intersect but not valid relation
  760. (Point(1, 1), []),
  761. # all points of result from tree are in common with box
  762. (box(3, 3, 6, 6), []),
  763. # all points of result from tree are in common with buffer
  764. (shapely.buffer(Point(3, 3), 1), []),
  765. # only one point of multipoint intersects but not valid relation
  766. (MultiPoint([[5, 7], [7, 7]]), []),
  767. ],
  768. )
  769. def test_query_crosses_points(tree, geometry, expected):
  770. assert_array_equal(tree.query(geometry, predicate="crosses"), expected)
  771. @pytest.mark.parametrize(
  772. "geometry,expected",
  773. [
  774. # point intersects first line but is completely in common with line
  775. (Point(0, 0), []),
  776. ([Point(0, 0)], [[], []]),
  777. # box overlaps envelope of first 2 lines, contains first and crosses second
  778. (box(0, 0, 1.5, 1.5), [1]),
  779. ([box(0, 0, 1.5, 1.5)], [[0], [1]]),
  780. # buffer intersects 2 lines
  781. (shapely.buffer(Point(3, 3), 0.5), [2, 3]),
  782. ([shapely.buffer(Point(3, 3), 0.5)], [[0, 0], [2, 3]]),
  783. # line crosses line
  784. (shapely.linestrings([(1, 0), (0, 1)]), [0]),
  785. ([shapely.linestrings([(1, 0), (0, 1)])], [[0], [0]]),
  786. # envelope of points overlaps lines but intersects none
  787. (MultiPoint([[5, 7], [7, 5]]), []),
  788. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  789. # only one point of multipoint intersects
  790. (MultiPoint([[5, 7], [7, 7], [7, 8]]), []),
  791. ([MultiPoint([[5, 7], [7, 7], [7, 8]])], [[], []]),
  792. ],
  793. )
  794. def test_query_crosses_lines(line_tree, geometry, expected):
  795. assert_array_equal(line_tree.query(geometry, predicate="crosses"), expected)
  796. @pytest.mark.parametrize(
  797. "geometry,expected",
  798. [
  799. # point within first polygon but not valid relation
  800. (Point(0, 0.5), []),
  801. ([Point(0, 0.5)], [[], []]),
  802. # box overlaps 2 polygons but not valid relation
  803. (box(0, 0, 1.5, 1.5), []),
  804. ([box(0, 0, 1.5, 1.5)], [[], []]),
  805. # buffer overlaps 3 polygons but not valid relation
  806. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), []),
  807. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[], []]),
  808. # only one point of multipoint within
  809. (MultiPoint([[5, 7], [7, 7], [7, 8]]), [7]),
  810. ([MultiPoint([[5, 7], [7, 7], [7, 8]])], [[0], [7]]),
  811. ],
  812. )
  813. def test_query_crosses_polygons(poly_tree, geometry, expected):
  814. assert_array_equal(poly_tree.query(geometry, predicate="crosses"), expected)
  815. ### predicate == 'touches'
  816. # See: https://postgis.net/docs/ST_Touches.html
  817. @pytest.mark.parametrize(
  818. "geometry,expected",
  819. [
  820. # points do not intersect
  821. (Point(0.5, 0.5), []),
  822. ([Point(0.5, 0.5)], [[], []]),
  823. # points intersect but not valid relation
  824. (Point(1, 1), []),
  825. ([Point(1, 1)], [[], []]),
  826. # box contains points but touches only those at edges
  827. (box(3, 3, 6, 6), [3, 6]),
  828. ([box(3, 3, 6, 6)], [[0, 0], [3, 6]]),
  829. # polygon completely contains point in tree
  830. (shapely.buffer(Point(3, 3), 1), []),
  831. ([shapely.buffer(Point(3, 3), 1)], [[], []]),
  832. # linestring intersects 2 points but touches only one
  833. (LineString([(-1, -1), (1, 1)]), [1]),
  834. ([LineString([(-1, -1), (1, 1)])], [[0], [1]]),
  835. # multipoints intersect but not valid relation
  836. (MultiPoint([[5, 5], [7, 7]]), []),
  837. ([MultiPoint([[5, 5], [7, 7]])], [[], []]),
  838. ],
  839. )
  840. def test_query_touches_points(tree, geometry, expected):
  841. assert_array_equal(tree.query(geometry, predicate="touches"), expected)
  842. @pytest.mark.parametrize(
  843. "geometry,expected",
  844. [
  845. # point intersects first line
  846. (Point(0, 0), [0]),
  847. ([Point(0, 0)], [[0], [0]]),
  848. # point is within line
  849. (Point(0.5, 0.5), []),
  850. ([Point(0.5, 0.5)], [[], []]),
  851. # point at shared vertex between 2 lines
  852. (Point(1, 1), [0, 1]),
  853. ([Point(1, 1)], [[0, 0], [0, 1]]),
  854. # box overlaps envelope of first 2 lines (touches edge of 1)
  855. (box(0, 0, 1, 1), [1]),
  856. ([box(0, 0, 1, 1)], [[0], [1]]),
  857. # buffer intersects 2 lines but does not touch edges of either
  858. (shapely.buffer(Point(3, 3), 0.5), []),
  859. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  860. # buffer intersects midpoint of line at tangent but there is a little overlap
  861. # due to precision issues
  862. (shapely.buffer(Point(2, 1), HALF_UNIT_DIAG + 1e-7), []),
  863. ([shapely.buffer(Point(2, 1), HALF_UNIT_DIAG + 1e-7)], [[], []]),
  864. # envelope of points overlaps lines but intersects none
  865. (MultiPoint([[5, 7], [7, 5]]), []),
  866. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  867. # only one point of multipoint intersects at vertex between lines
  868. (MultiPoint([[5, 7], [7, 7], [7, 8]]), [6, 7]),
  869. ([MultiPoint([[5, 7], [7, 7], [7, 8]])], [[0, 0], [6, 7]]),
  870. ],
  871. )
  872. def test_query_touches_lines(line_tree, geometry, expected):
  873. assert_array_equal(line_tree.query(geometry, predicate="touches"), expected)
  874. @pytest.mark.parametrize(
  875. "geometry,expected",
  876. [
  877. # point within first polygon
  878. (Point(0, 0.5), []),
  879. ([Point(0, 0.5)], [[], []]),
  880. # point is at edge of first polygon
  881. (Point(HALF_UNIT_DIAG + EPS, 0), [0]),
  882. ([Point(HALF_UNIT_DIAG + EPS, 0)], [[0], [0]]),
  883. # box overlaps envelope of 2 polygons does not touch any at edge
  884. (box(0, 0, 1, 1), []),
  885. ([box(0, 0, 1, 1)], [[], []]),
  886. # box overlaps 2 polygons and touches edge of first
  887. (box(HALF_UNIT_DIAG + EPS, 0, 2, 2), [0]),
  888. ([box(HALF_UNIT_DIAG + EPS, 0, 2, 2)], [[0], [0]]),
  889. # buffer overlaps 3 polygons but does not touch any at edge
  890. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG + EPS), []),
  891. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG + EPS)], [[], []]),
  892. # only one point of multipoint within polygon but does not touch
  893. (MultiPoint([[0, 0], [7, 7], [7, 8]]), []),
  894. ([MultiPoint([[0, 0], [7, 7], [7, 8]])], [[], []]),
  895. ],
  896. )
  897. def test_query_touches_polygons(poly_tree, geometry, expected):
  898. assert_array_equal(poly_tree.query(geometry, predicate="touches"), expected)
  899. ### predicate == 'covers'
  900. @pytest.mark.parametrize(
  901. "geometry,expected",
  902. [
  903. # points do not intersect
  904. (Point(0.5, 0.5), []),
  905. ([Point(0.5, 0.5)], [[], []]),
  906. # points intersect and thus no point is outside the other
  907. (Point(1, 1), [1]),
  908. ([Point(1, 1)], [[0], [1]]),
  909. # box covers any points that intersect or are within
  910. (box(3, 3, 6, 6), [3, 4, 5, 6]),
  911. ([box(3, 3, 6, 6)], [[0, 0, 0, 0], [3, 4, 5, 6]]),
  912. # envelope of buffer covers more points than are covered by buffer
  913. # due to diagonal distance
  914. (shapely.buffer(Point(3, 3), 1), [3]),
  915. ([shapely.buffer(Point(3, 3), 1)], [[0], [3]]),
  916. # envelope of buffer with 1/2 distance between points should intersect
  917. # same points as envelope
  918. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [2, 3, 4]),
  919. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[0, 0, 0], [2, 3, 4]]),
  920. # multipoints intersect and thus no point is outside the other
  921. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  922. ([MultiPoint([[5, 5], [7, 7]])], [[0, 0], [5, 7]]),
  923. # envelope of points contains points, but points do not intersect
  924. (MultiPoint([[5, 7], [7, 5]]), []),
  925. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  926. # only one point of multipoint intersects
  927. (MultiPoint([[5, 7], [7, 7]]), [7]),
  928. ([MultiPoint([[5, 7], [7, 7]])], [[0], [7]]),
  929. ],
  930. )
  931. def test_query_covers_points(tree, geometry, expected):
  932. assert_array_equal(tree.query(geometry, predicate="covers"), expected)
  933. @pytest.mark.parametrize(
  934. "geometry,expected",
  935. [
  936. # point does not cover any lines (not valid relation)
  937. (Point(0, 0), []),
  938. ([Point(0, 0)], [[], []]),
  939. # box covers first line (intersects another does not contain it)
  940. (box(0, 0, 1.5, 1.5), [0]),
  941. ([box(0, 0, 1.5, 1.5)], [[0], [0]]),
  942. # box completely covers 2 lines (touches edges of 2 others)
  943. (box(1, 1, 3, 3), [1, 2]),
  944. ([box(1, 1, 3, 3)], [[0, 0], [1, 2]]),
  945. # buffer intersects 2 lines but does not completely cover either
  946. (shapely.buffer(Point(3, 3), 0.5), []),
  947. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  948. # envelope of points overlaps lines but intersects none
  949. (MultiPoint([[5, 7], [7, 5]]), []),
  950. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  951. # only one point of multipoint intersects a line, but does not completely cover
  952. # it
  953. (MultiPoint([[5, 7], [7, 7]]), []),
  954. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  955. # both points intersect but do not cover any lines (not valid relation)
  956. (MultiPoint([[5, 5], [6, 6]]), []),
  957. ([MultiPoint([[5, 5], [6, 6]])], [[], []]),
  958. ],
  959. )
  960. def test_query_covers_lines(line_tree, geometry, expected):
  961. assert_array_equal(line_tree.query(geometry, predicate="covers"), expected)
  962. @pytest.mark.parametrize(
  963. "geometry,expected",
  964. [
  965. # point does not cover any polygons (not valid relation)
  966. (Point(0, 0), []),
  967. ([Point(0, 0)], [[], []]),
  968. # box overlaps envelope of 2 polygons but does not completely cover either
  969. (box(0, 0, 1, 1), []),
  970. ([box(0, 0, 1, 1)], [[], []]),
  971. # larger box intersects 3 polygons but covers only one
  972. (box(0, 0, 2, 2), [1]),
  973. ([box(0, 0, 2, 2)], [[0], [1]]),
  974. # buffer overlaps 3 polygons but does not completely cover any
  975. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), []),
  976. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[], []]),
  977. # larger buffer overlaps 6 polygons (touches midpoints) but covers only one
  978. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [3]),
  979. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[0], [3]]),
  980. # envelope of points overlaps polygons, but points do not intersect
  981. # (not valid relation)
  982. (MultiPoint([[5, 7], [7, 5]]), []),
  983. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  984. ],
  985. )
  986. def test_query_covers_polygons(poly_tree, geometry, expected):
  987. assert_array_equal(poly_tree.query(geometry, predicate="covers"), expected)
  988. ### predicate == 'covered_by'
  989. @pytest.mark.parametrize(
  990. "geometry,expected",
  991. [
  992. # points do not intersect
  993. (Point(0.5, 0.5), []),
  994. ([Point(0.5, 0.5)], [[], []]),
  995. # points intersect
  996. (Point(1, 1), [1]),
  997. ([Point(1, 1)], [[0], [1]]),
  998. # box not covered by points
  999. (box(3, 3, 6, 6), []),
  1000. ([box(3, 3, 6, 6)], [[], []]),
  1001. # envelope of buffer not covered by points
  1002. (shapely.buffer(Point(3, 3), 1), []),
  1003. ([shapely.buffer(Point(3, 3), 1)], [[], []]),
  1004. # multipoints intersect but are not covered by points in tree
  1005. (MultiPoint([[5, 5], [7, 7]]), []),
  1006. ([MultiPoint([[5, 5], [7, 7]])], [[], []]),
  1007. # only one point of multipoint intersects, but multipoints are not
  1008. # covered by any points in tree
  1009. (MultiPoint([[5, 7], [7, 7]]), []),
  1010. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  1011. # envelope of points overlaps points, but points do not intersect
  1012. (MultiPoint([[5, 7], [7, 5]]), []),
  1013. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  1014. ],
  1015. )
  1016. def test_query_covered_by_points(tree, geometry, expected):
  1017. assert_array_equal(tree.query(geometry, predicate="covered_by"), expected)
  1018. @pytest.mark.parametrize(
  1019. "geometry,expected",
  1020. [
  1021. # endpoint is covered by first line
  1022. (Point(0, 0), [0]),
  1023. ([Point(0, 0)], [[0], [0]]),
  1024. # point covered by first line
  1025. (Point(0.5, 0.5), [0]),
  1026. ([Point(0.5, 0.5)], [[0], [0]]),
  1027. # point within envelope of first line but does not intersect
  1028. (Point(0, 0.5), []),
  1029. ([Point(0, 0.5)], [[], []]),
  1030. # point at shared vertex between 2 lines and is covered by both
  1031. (Point(1, 1), [0, 1]),
  1032. ([Point(1, 1)], [[0, 0], [0, 1]]),
  1033. # line intersects 3 lines, but is covered by only one
  1034. (shapely.linestrings([[1, 1], [2, 2]]), [1]),
  1035. ([shapely.linestrings([[1, 1], [2, 2]])], [[0], [1]]),
  1036. # line intersects 2 lines, but is covered by neither
  1037. (shapely.linestrings([[1.5, 1.5], [2.5, 2.5]]), []),
  1038. ([shapely.linestrings([[1.5, 1.5], [2.5, 2.5]])], [[], []]),
  1039. # box not covered by line (not valid geometric relation)
  1040. (box(0, 0, 1, 1), []),
  1041. ([box(0, 0, 1, 1)], [[], []]),
  1042. # buffer intersects 2 lines but not within either (not valid geometric relation)
  1043. (shapely.buffer(Point(3, 3), 0.5), []),
  1044. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  1045. # envelope of points overlaps lines but intersects none
  1046. (MultiPoint([[5, 7], [7, 5]]), []),
  1047. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  1048. # only one point of multipoint intersects, but both are not covered by line
  1049. (MultiPoint([[5, 7], [7, 7]]), []),
  1050. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  1051. # both points are covered by a line
  1052. (MultiPoint([[6.5, 6.5], [7, 7]]), [6]),
  1053. ([MultiPoint([[6.5, 6.5], [7, 7]])], [[0], [6]]),
  1054. ],
  1055. )
  1056. def test_query_covered_by_lines(line_tree, geometry, expected):
  1057. assert_array_equal(line_tree.query(geometry, predicate="covered_by"), expected)
  1058. @pytest.mark.parametrize(
  1059. "geometry,expected",
  1060. [
  1061. # point covered by polygon
  1062. (Point(0, 0.5), [0]),
  1063. ([Point(0, 0.5)], [[0], [0]]),
  1064. (Point(0.5, 0), [0]),
  1065. ([Point(0.5, 0)], [[0], [0]]),
  1066. (Point(1, 1), [1]),
  1067. ([Point(1, 1)], [[0], [1]]),
  1068. # midpoint between two polygons is covered by both
  1069. (Point(0.5, 0.5), [0, 1]),
  1070. ([Point(0.5, 0.5)], [[0, 0], [0, 1]]),
  1071. # line intersects multiple polygons but is not covered by any
  1072. (shapely.linestrings([[0, 0], [2, 2]]), []),
  1073. ([shapely.linestrings([[0, 0], [2, 2]])], [[], []]),
  1074. # line intersects multiple polygons but is covered by only one
  1075. (shapely.linestrings([[1.5, 1.5], [2.5, 2.5]]), [2]),
  1076. ([shapely.linestrings([[1.5, 1.5], [2.5, 2.5]])], [[0], [2]]),
  1077. # box overlaps envelope of 2 polygons but not covered by either
  1078. (box(0, 0, 1, 1), []),
  1079. ([box(0, 0, 1, 1)], [[], []]),
  1080. # box covered by polygon
  1081. (box(0, 0, 0.5, 0.5), [0]),
  1082. ([box(0, 0, 0.5, 0.5)], [[0], [0]]),
  1083. # larger box intersects 3 polygons but not covered by any
  1084. (box(0, 0, 1.5, 1.5), []),
  1085. ([box(0, 0, 1.5, 1.5)], [[], []]),
  1086. # buffer intersects 3 polygons but only within one
  1087. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [3]),
  1088. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[0], [3]]),
  1089. # larger buffer overlaps 6 polygons (touches midpoints) but within none
  1090. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), []),
  1091. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[], []]),
  1092. # envelope of points overlaps polygons, but points do not intersect
  1093. (MultiPoint([[5, 7], [7, 5]]), []),
  1094. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  1095. # only one point of multipoint within polygon
  1096. (MultiPoint([[5, 7], [7, 7]]), []),
  1097. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  1098. # both points in multipoint within polygon
  1099. (MultiPoint([[5.25, 5.5], [5.25, 5.0]]), [5]),
  1100. ([MultiPoint([[5.25, 5.5], [5.25, 5.0]])], [[0], [5]]),
  1101. ],
  1102. )
  1103. def test_query_covered_by_polygons(poly_tree, geometry, expected):
  1104. assert_array_equal(poly_tree.query(geometry, predicate="covered_by"), expected)
  1105. ### predicate == 'contains_properly'
  1106. @pytest.mark.parametrize(
  1107. "geometry,expected",
  1108. [
  1109. # points do not intersect
  1110. (Point(0.5, 0.5), []),
  1111. ([Point(0.5, 0.5)], [[], []]),
  1112. # points intersect
  1113. (Point(1, 1), [1]),
  1114. ([Point(1, 1)], [[0], [1]]),
  1115. # line contains every point that is not on its first or last coordinate
  1116. # these are on the "exterior" of the line
  1117. (shapely.linestrings([[0, 0], [2, 2]]), [1]),
  1118. ([shapely.linestrings([[0, 0], [2, 2]])], [[0], [1]]),
  1119. # slightly longer line contains multiple points
  1120. (shapely.linestrings([[0.5, 0.5], [2.5, 2.5]]), [1, 2]),
  1121. ([shapely.linestrings([[0.5, 0.5], [2.5, 2.5]])], [[0, 0], [1, 2]]),
  1122. # line intersects and contains one point
  1123. (shapely.linestrings([[0, 2], [2, 0]]), [1]),
  1124. ([shapely.linestrings([[0, 2], [2, 0]])], [[0], [1]]),
  1125. # box contains points (2 are at edges and not contained)
  1126. (box(3, 3, 6, 6), [4, 5]),
  1127. ([box(3, 3, 6, 6)], [[0, 0], [4, 5]]),
  1128. # envelope of buffer contains more points than within buffer
  1129. # due to diagonal distance
  1130. (shapely.buffer(Point(3, 3), 1), [3]),
  1131. ([shapely.buffer(Point(3, 3), 1)], [[0], [3]]),
  1132. # envelope of buffer with 1/2 distance between points should intersect
  1133. # same points as envelope
  1134. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [2, 3, 4]),
  1135. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[0, 0, 0], [2, 3, 4]]),
  1136. # multipoints intersect
  1137. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  1138. ([MultiPoint([[5, 5], [7, 7]])], [[0, 0], [5, 7]]),
  1139. # envelope of points contains points, but points do not intersect
  1140. (MultiPoint([[5, 7], [7, 5]]), []),
  1141. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  1142. # only one point of multipoint intersects
  1143. (MultiPoint([[5, 7], [7, 7]]), [7]),
  1144. ([MultiPoint([[5, 7], [7, 7]])], [[0], [7]]),
  1145. ],
  1146. )
  1147. def test_query_contains_properly_points(tree, geometry, expected):
  1148. assert_array_equal(tree.query(geometry, predicate="contains_properly"), expected)
  1149. @pytest.mark.parametrize(
  1150. "geometry,expected",
  1151. [
  1152. # None of the following conditions satisfy the relation for linestrings
  1153. # because they have no interior:
  1154. # "a contains b if no points of b lie in the exterior of a, and at least one
  1155. # point of the interior of b lies in the interior of a"
  1156. (Point(0, 0), []),
  1157. ([Point(0, 0)], [[], []]),
  1158. (shapely.linestrings([[0, 0], [1, 1]]), []),
  1159. ([shapely.linestrings([[0, 0], [1, 1]])], [[], []]),
  1160. (shapely.linestrings([[0, 0], [2, 2]]), []),
  1161. ([shapely.linestrings([[0, 0], [2, 2]])], [[], []]),
  1162. (shapely.linestrings([[0, 2], [2, 0]]), []),
  1163. ([shapely.linestrings([[0, 2], [2, 0]])], [[], []]),
  1164. (MultiPoint([[5, 7], [7, 5]]), []),
  1165. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  1166. (MultiPoint([[5, 7], [7, 7]]), []),
  1167. ([MultiPoint([[5, 7], [7, 7]])], [[], []]),
  1168. (MultiPoint([[5, 5], [6, 6]]), []),
  1169. ([MultiPoint([[5, 5], [6, 6]])], [[], []]),
  1170. (box(0, 0, 1, 1), []),
  1171. ([box(0, 0, 1, 1)], [[], []]),
  1172. (box(0, 0, 2, 2), []),
  1173. ([box(0, 0, 2, 2)], [[], []]),
  1174. (shapely.buffer(Point(3, 3), 0.5), []),
  1175. ([shapely.buffer(Point(3, 3), 0.5)], [[], []]),
  1176. ],
  1177. )
  1178. def test_query_contains_properly_lines(line_tree, geometry, expected):
  1179. assert_array_equal(
  1180. line_tree.query(geometry, predicate="contains_properly"), expected
  1181. )
  1182. @pytest.mark.parametrize(
  1183. "geometry,expected",
  1184. [
  1185. # point does not contain any polygons (not valid relation)
  1186. (Point(0, 0), []),
  1187. ([Point(0, 0)], [[], []]),
  1188. # line intersects multiple polygons but does not contain any (not valid
  1189. # relation)
  1190. (shapely.linestrings([[0, 0], [2, 2]]), []),
  1191. ([shapely.linestrings([[0, 0], [2, 2]])], [[], []]),
  1192. # box overlaps envelope of 2 polygons but contains neither
  1193. (box(0, 0, 1, 1), []),
  1194. ([box(0, 0, 1, 1)], [[], []]),
  1195. # larger box intersects 3 polygons but contains only one
  1196. (box(0, 0, 2, 2), [1]),
  1197. ([box(0, 0, 2, 2)], [[0], [1]]),
  1198. # buffer overlaps 3 polygons but contains none
  1199. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), []),
  1200. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[], []]),
  1201. # larger buffer overlaps 6 polygons (touches midpoints) but contains one
  1202. (shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG), [3]),
  1203. ([shapely.buffer(Point(3, 3), 3 * HALF_UNIT_DIAG)], [[0], [3]]),
  1204. # envelope of points overlaps polygons, but points do not intersect
  1205. # (not valid relation)
  1206. (MultiPoint([[5, 7], [7, 5]]), []),
  1207. ([MultiPoint([[5, 7], [7, 5]])], [[], []]),
  1208. ],
  1209. )
  1210. def test_query_contains_properly_polygons(poly_tree, geometry, expected):
  1211. assert_array_equal(
  1212. poly_tree.query(geometry, predicate="contains_properly"), expected
  1213. )
  1214. ### predicate = 'dwithin'
  1215. @pytest.mark.skipif(geos_version >= (3, 10, 0), reason="GEOS >= 3.10")
  1216. @pytest.mark.parametrize(
  1217. "geometry", [Point(0, 0), [Point(0, 0)], None, [None], empty, [empty]]
  1218. )
  1219. def test_query_dwithin_geos_version(tree, geometry):
  1220. with pytest.raises(UnsupportedGEOSVersionError, match="requires GEOS >= 3.10"):
  1221. tree.query(geometry, predicate="dwithin", distance=1)
  1222. @pytest.mark.skipif(geos_version < (3, 10, 0), reason="GEOS < 3.10")
  1223. @pytest.mark.parametrize(
  1224. "geometry,distance,match",
  1225. [
  1226. (Point(0, 0), None, "distance parameter must be provided"),
  1227. ([Point(0, 0)], None, "distance parameter must be provided"),
  1228. (Point(0, 0), "foo", "could not convert string to float"),
  1229. ([Point(0, 0)], "foo", "could not convert string to float"),
  1230. ([Point(0, 0)], ["foo"], "could not convert string to float"),
  1231. (Point(0, 0), [0, 1], "Could not broadcast distance to match geometry"),
  1232. ([Point(0, 0)], [0, 1], "Could not broadcast distance to match geometry"),
  1233. (Point(0, 0), [[1.0]], "should be one dimensional"),
  1234. ([Point(0, 0)], [[1.0]], "should be one dimensional"),
  1235. ],
  1236. )
  1237. def test_query_dwithin_invalid_distance(tree, geometry, distance, match):
  1238. with pytest.raises(ValueError, match=match):
  1239. tree.query(geometry, predicate="dwithin", distance=distance)
  1240. @pytest.mark.skipif(geos_version < (3, 10, 0), reason="GEOS < 3.10")
  1241. @pytest.mark.parametrize(
  1242. "geometry,distance,expected",
  1243. [
  1244. (None, 1.0, []),
  1245. ([None], 1.0, [[], []]),
  1246. (Point(0.25, 0.25), 0, []),
  1247. ([Point(0.25, 0.25)], 0, [[], []]),
  1248. (Point(0.25, 0.25), -1, []),
  1249. ([Point(0.25, 0.25)], -1, [[], []]),
  1250. (Point(0.25, 0.25), np.nan, []),
  1251. ([Point(0.25, 0.25)], np.nan, [[], []]),
  1252. (Point(), 1, []),
  1253. ([Point()], 1, [[], []]),
  1254. (Point(0.25, 0.25), 0.5, [0]),
  1255. ([Point(0.25, 0.25)], 0.5, [[0], [0]]),
  1256. (Point(0.25, 0.25), 2.5, [0, 1, 2]),
  1257. ([Point(0.25, 0.25)], 2.5, [[0, 0, 0], [0, 1, 2]]),
  1258. (Point(3, 3), 1.5, [2, 3, 4]),
  1259. ([Point(3, 3)], 1.5, [[0, 0, 0], [2, 3, 4]]),
  1260. # 2 equidistant points in tree
  1261. (Point(0.5, 0.5), 0.75, [0, 1]),
  1262. ([Point(0.5, 0.5)], 0.75, [[0, 0], [0, 1]]),
  1263. (
  1264. [None, Point(0.5, 0.5)],
  1265. 0.75,
  1266. [
  1267. [
  1268. 1,
  1269. 1,
  1270. ],
  1271. [0, 1],
  1272. ],
  1273. ),
  1274. (
  1275. [Point(0.5, 0.5), Point(0.25, 0.25)],
  1276. 0.75,
  1277. [[0, 0, 1], [0, 1, 0]],
  1278. ),
  1279. (
  1280. [Point(0, 0.2), Point(1.75, 1.75)],
  1281. [0.25, 2],
  1282. [[0, 1, 1, 1], [0, 1, 2, 3]],
  1283. ),
  1284. # all points intersect box
  1285. (box(0, 0, 3, 3), 0, [0, 1, 2, 3]),
  1286. ([box(0, 0, 3, 3)], 0, [[0, 0, 0, 0], [0, 1, 2, 3]]),
  1287. (box(0, 0, 3, 3), 0.25, [0, 1, 2, 3]),
  1288. ([box(0, 0, 3, 3)], 0.25, [[0, 0, 0, 0], [0, 1, 2, 3]]),
  1289. # intersecting and nearby points
  1290. (box(1, 1, 2, 2), 1.5, [0, 1, 2, 3]),
  1291. ([box(1, 1, 2, 2)], 1.5, [[0, 0, 0, 0], [0, 1, 2, 3]]),
  1292. # # return nearest point in tree for each point in multipoint
  1293. (MultiPoint([[0.25, 0.25], [1.5, 1.5]]), 0.75, [0, 1, 2]),
  1294. ([MultiPoint([[0.25, 0.25], [1.5, 1.5]])], 0.75, [[0, 0, 0], [0, 1, 2]]),
  1295. # 2 equidistant points per point in multipoint
  1296. (
  1297. MultiPoint([[0.5, 0.5], [3.5, 3.5]]),
  1298. 0.75,
  1299. [0, 1, 3, 4],
  1300. ),
  1301. (
  1302. [MultiPoint([[0.5, 0.5], [3.5, 3.5]])],
  1303. 0.75,
  1304. [[0, 0, 0, 0], [0, 1, 3, 4]],
  1305. ),
  1306. ],
  1307. )
  1308. def test_query_dwithin_points(tree, geometry, distance, expected):
  1309. assert_array_equal(
  1310. tree.query(geometry, predicate="dwithin", distance=distance), expected
  1311. )
  1312. @pytest.mark.skipif(geos_version < (3, 10, 0), reason="GEOS < 3.10")
  1313. @pytest.mark.parametrize(
  1314. "geometry,distance,expected",
  1315. [
  1316. (None, 1.0, []),
  1317. ([None], 1.0, [[], []]),
  1318. (Point(0.5, 0.5), 0, [0]),
  1319. ([Point(0.5, 0.5)], 0, [[0], [0]]),
  1320. (Point(0.5, 0.5), 1.0, [0, 1]),
  1321. ([Point(0.5, 0.5)], 1.0, [[0, 0], [0, 1]]),
  1322. (Point(2, 2), 0.5, [1, 2]),
  1323. ([Point(2, 2)], 0.5, [[0, 0], [1, 2]]),
  1324. (box(0, 0, 1, 1), 0.5, [0, 1]),
  1325. ([box(0, 0, 1, 1)], 0.5, [[0, 0], [0, 1]]),
  1326. (box(0.5, 0.5, 1.5, 1.5), 0.5, [0, 1]),
  1327. ([box(0.5, 0.5, 1.5, 1.5)], 0.5, [[0, 0], [0, 1]]),
  1328. # multipoints at endpoints of 2 lines each
  1329. (MultiPoint([[5, 5], [7, 7]]), 0.5, [4, 5, 6, 7]),
  1330. ([MultiPoint([[5, 5], [7, 7]])], 0.5, [[0, 0, 0, 0], [4, 5, 6, 7]]),
  1331. # multipoints are equidistant from 2 lines
  1332. (MultiPoint([[5, 7], [7, 5]]), 1.5, [5, 6]),
  1333. ([MultiPoint([[5, 7], [7, 5]])], 1.5, [[0, 0], [5, 6]]),
  1334. ],
  1335. )
  1336. def test_query_dwithin_lines(line_tree, geometry, distance, expected):
  1337. assert_array_equal(
  1338. line_tree.query(geometry, predicate="dwithin", distance=distance),
  1339. expected,
  1340. )
  1341. @pytest.mark.skipif(geos_version < (3, 10, 0), reason="GEOS < 3.10")
  1342. @pytest.mark.parametrize(
  1343. "geometry,distance,expected",
  1344. [
  1345. (Point(0, 0), 0, [0]),
  1346. ([Point(0, 0)], 0, [[0], [0]]),
  1347. (Point(0, 0), 0.5, [0]),
  1348. ([Point(0, 0)], 0.5, [[0], [0]]),
  1349. (Point(0, 0), 1.5, [0, 1]),
  1350. ([Point(0, 0)], 1.5, [[0, 0], [0, 1]]),
  1351. (Point(0.5, 0.5), 1, [0, 1]),
  1352. ([Point(0.5, 0.5)], 1, [[0, 0], [0, 1]]),
  1353. (Point(0.5, 0.5), 0.5, [0, 1]),
  1354. ([Point(0.5, 0.5)], 0.5, [[0, 0], [0, 1]]),
  1355. (box(0, 0, 1, 1), 0, [0, 1]),
  1356. ([box(0, 0, 1, 1)], 0, [[0, 0], [0, 1]]),
  1357. (box(0, 0, 1, 1), 2, [0, 1, 2]),
  1358. ([box(0, 0, 1, 1)], 2, [[0, 0, 0], [0, 1, 2]]),
  1359. (MultiPoint([[5, 5], [7, 7]]), 0.5, [5, 7]),
  1360. ([MultiPoint([[5, 5], [7, 7]])], 0.5, [[0, 0], [5, 7]]),
  1361. (
  1362. MultiPoint([[5, 5], [7, 7]]),
  1363. 2.5,
  1364. [3, 4, 5, 6, 7, 8, 9],
  1365. ),
  1366. (
  1367. [MultiPoint([[5, 5], [7, 7]])],
  1368. 2.5,
  1369. [[0, 0, 0, 0, 0, 0, 0], [3, 4, 5, 6, 7, 8, 9]],
  1370. ),
  1371. ],
  1372. )
  1373. def test_query_dwithin_polygons(poly_tree, geometry, distance, expected):
  1374. assert_array_equal(
  1375. poly_tree.query(geometry, predicate="dwithin", distance=distance),
  1376. expected,
  1377. )
  1378. ### STRtree nearest
  1379. def test_nearest_empty_tree():
  1380. tree = STRtree([])
  1381. assert tree.nearest(point) is None
  1382. @pytest.mark.parametrize("geometry", ["I am not a geometry"])
  1383. def test_nearest_invalid_geom(tree, geometry):
  1384. with pytest.raises(TypeError):
  1385. tree.nearest(geometry)
  1386. @pytest.mark.parametrize("geometry", [None, [None], [Point(1, 1), None]])
  1387. def test_nearest_none(tree, geometry):
  1388. with pytest.raises(ValueError):
  1389. tree.nearest(geometry)
  1390. @pytest.mark.parametrize(
  1391. "geometry", [empty_point, [empty_point], [Point(1, 1), empty_point]]
  1392. )
  1393. def test_nearest_empty(tree, geometry):
  1394. with pytest.raises(ValueError):
  1395. tree.nearest(geometry)
  1396. @pytest.mark.parametrize(
  1397. "geometry,expected",
  1398. [
  1399. (Point(0.25, 0.25), 0),
  1400. (Point(0.75, 0.75), 1),
  1401. (Point(1, 1), 1),
  1402. ([Point(1, 1), Point(0, 0)], [1, 0]),
  1403. ([Point(1, 1), Point(0.25, 1)], [1, 1]),
  1404. ([Point(-10, -10), Point(100, 100)], [0, 9]),
  1405. (box(0.5, 0.5, 0.75, 0.75), 1),
  1406. (shapely.buffer(Point(2.5, 2.5), HALF_UNIT_DIAG), 2),
  1407. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), 3),
  1408. (MultiPoint([[5.5, 5], [7, 7]]), 7),
  1409. (MultiPoint([[5, 7], [7, 5]]), 6),
  1410. ],
  1411. )
  1412. def test_nearest_points(tree, geometry, expected):
  1413. assert_array_equal(tree.nearest(geometry), expected)
  1414. @pytest.mark.parametrize(
  1415. "geometry,expected",
  1416. [
  1417. # 2 equidistant points in tree
  1418. (Point(0.5, 0.5), [0, 1]),
  1419. # multiple points in box
  1420. (box(0, 0, 3, 3), [0, 1, 2, 3]),
  1421. # return nearest point in tree for each point in multipoint
  1422. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  1423. ],
  1424. )
  1425. def test_nearest_points_equidistant(tree, geometry, expected):
  1426. # results are returned in order they are traversed when searching the tree,
  1427. # which can vary between GEOS versions, so we test that one of the valid
  1428. # results is present
  1429. result = tree.nearest(geometry)
  1430. assert result in expected
  1431. @pytest.mark.parametrize(
  1432. "geometry,expected",
  1433. [
  1434. (Point(0.5, 0.5), 0),
  1435. (Point(1.5, 0.5), 0),
  1436. (shapely.box(0.5, 1.5, 1, 2), 1),
  1437. (shapely.linestrings([[0, 0.5], [1, 2.5]]), 0),
  1438. ],
  1439. )
  1440. def test_nearest_lines(line_tree, geometry, expected):
  1441. assert_array_equal(line_tree.nearest(geometry), expected)
  1442. @pytest.mark.parametrize(
  1443. "geometry,expected",
  1444. [
  1445. # at junction between 2 lines
  1446. (Point(2, 2), [1, 2]),
  1447. # contains one line, intersects with another
  1448. (box(0, 0, 1, 1), [0, 1]),
  1449. # overlaps 2 lines
  1450. (box(0.5, 0.5, 1.5, 1.5), [0, 1]),
  1451. # box overlaps 2 lines and intersects endpoints of 2 more
  1452. (box(3, 3, 5, 5), [2, 3, 4, 5]),
  1453. (shapely.buffer(Point(2.5, 2.5), HALF_UNIT_DIAG), [1, 2]),
  1454. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [2, 3]),
  1455. # multipoints at endpoints of 2 lines each
  1456. (MultiPoint([[5, 5], [7, 7]]), [4, 5, 6, 7]),
  1457. # second point in multipoint at endpoints of 2 lines
  1458. (MultiPoint([[5.5, 5], [7, 7]]), [6, 7]),
  1459. # multipoints are equidistant from 2 lines
  1460. (MultiPoint([[5, 7], [7, 5]]), [5, 6]),
  1461. ],
  1462. )
  1463. def test_nearest_lines_equidistant(line_tree, geometry, expected):
  1464. # results are returned in order they are traversed when searching the tree,
  1465. # which can vary between GEOS versions, so we test that one of the valid
  1466. # results is present
  1467. result = line_tree.nearest(geometry)
  1468. assert result in expected
  1469. @pytest.mark.parametrize(
  1470. "geometry,expected",
  1471. [
  1472. (Point(0, 0), 0),
  1473. (Point(2, 2), 2),
  1474. (shapely.box(0, 5, 1, 6), 3),
  1475. (MultiPoint([[5, 7], [7, 5]]), 6),
  1476. ],
  1477. )
  1478. def test_nearest_polygons(poly_tree, geometry, expected):
  1479. assert_array_equal(poly_tree.nearest(geometry), expected)
  1480. @pytest.mark.parametrize(
  1481. "geometry,expected",
  1482. [
  1483. # 2 polygons in tree overlap point
  1484. (Point(0.5, 0.5), [0, 1]),
  1485. # box overlaps multiple polygons
  1486. (box(0, 0, 1, 1), [0, 1]),
  1487. (box(0.5, 0.5, 1.5, 1.5), [0, 1, 2]),
  1488. (box(3, 3, 5, 5), [3, 4, 5]),
  1489. (shapely.buffer(Point(2.5, 2.5), HALF_UNIT_DIAG), [2, 3]),
  1490. # completely overlaps one polygon, touches 2 others
  1491. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [2, 3, 4]),
  1492. # each point in multi point intersects a polygon in tree
  1493. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  1494. (MultiPoint([[5.5, 5], [7, 7]]), [5, 7]),
  1495. ],
  1496. )
  1497. def test_nearest_polygons_equidistant(poly_tree, geometry, expected):
  1498. # results are returned in order they are traversed when searching the tree,
  1499. # which can vary between GEOS versions, so we test that one of the valid
  1500. # results is present
  1501. result = poly_tree.nearest(geometry)
  1502. assert result in expected
  1503. def test_query_nearest_empty_tree():
  1504. tree = STRtree([])
  1505. assert_array_equal(tree.query_nearest(point), [])
  1506. assert_array_equal(tree.query_nearest([point]), [[], []])
  1507. @pytest.mark.parametrize("geometry", ["I am not a geometry", ["still not a geometry"]])
  1508. def test_query_nearest_invalid_geom(tree, geometry):
  1509. with pytest.raises(TypeError):
  1510. tree.query_nearest(geometry)
  1511. @pytest.mark.parametrize(
  1512. "geometry,return_distance,expected",
  1513. [
  1514. (None, False, []),
  1515. ([None], False, [[], []]),
  1516. (None, True, ([], [])),
  1517. ([None], True, ([[], []], [])),
  1518. ],
  1519. )
  1520. def test_query_nearest_none(tree, geometry, return_distance, expected):
  1521. if return_distance:
  1522. index, distance = tree.query_nearest(geometry, return_distance=True)
  1523. assert_array_equal(index, expected[0])
  1524. assert_array_equal(distance, expected[1])
  1525. else:
  1526. assert_array_equal(tree.query_nearest(geometry), expected)
  1527. @pytest.mark.parametrize(
  1528. "geometry,expected",
  1529. [(empty, []), ([empty], [[], []]), ([empty, point], [[1, 1], [2, 3]])],
  1530. )
  1531. def test_query_nearest_empty_geom(tree, geometry, expected):
  1532. assert_array_equal(tree.query_nearest(geometry), expected)
  1533. @pytest.mark.parametrize(
  1534. "geometry,expected",
  1535. [
  1536. (Point(0.25, 0.25), [0]),
  1537. ([Point(0.25, 0.25)], [[0], [0]]),
  1538. (Point(0.75, 0.75), [1]),
  1539. ([Point(0.75, 0.75)], [[0], [1]]),
  1540. (Point(1, 1), [1]),
  1541. ([Point(1, 1)], [[0], [1]]),
  1542. # 2 equidistant points in tree
  1543. (Point(0.5, 0.5), [0, 1]),
  1544. ([Point(0.5, 0.5)], [[0, 0], [0, 1]]),
  1545. ([Point(1, 1), Point(0, 0)], [[0, 1], [1, 0]]),
  1546. ([Point(1, 1), Point(0.25, 1)], [[0, 1], [1, 1]]),
  1547. ([Point(-10, -10), Point(100, 100)], [[0, 1], [0, 9]]),
  1548. (box(0.5, 0.5, 0.75, 0.75), [1]),
  1549. ([box(0.5, 0.5, 0.75, 0.75)], [[0], [1]]),
  1550. # multiple points in box
  1551. (box(0, 0, 3, 3), [0, 1, 2, 3]),
  1552. ([box(0, 0, 3, 3)], [[0, 0, 0, 0], [0, 1, 2, 3]]),
  1553. (shapely.buffer(Point(2.5, 2.5), 1), [2, 3]),
  1554. ([shapely.buffer(Point(2.5, 2.5), 1)], [[0, 0], [2, 3]]),
  1555. (shapely.buffer(Point(3, 3), 0.5), [3]),
  1556. ([shapely.buffer(Point(3, 3), 0.5)], [[0], [3]]),
  1557. (MultiPoint([[5.5, 5], [7, 7]]), [7]),
  1558. ([MultiPoint([[5.5, 5], [7, 7]])], [[0], [7]]),
  1559. (MultiPoint([[5, 7], [7, 5]]), [6]),
  1560. ([MultiPoint([[5, 7], [7, 5]])], [[0], [6]]),
  1561. # return nearest point in tree for each point in multipoint
  1562. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  1563. ([MultiPoint([[5, 5], [7, 7]])], [[0, 0], [5, 7]]),
  1564. # 2 equidistant points per point in multipoint
  1565. (MultiPoint([[0.5, 0.5], [3.5, 3.5]]), [0, 1, 3, 4]),
  1566. ([MultiPoint([[0.5, 0.5], [3.5, 3.5]])], [[0, 0, 0, 0], [0, 1, 3, 4]]),
  1567. ],
  1568. )
  1569. def test_query_nearest_points(tree, geometry, expected):
  1570. assert_array_equal(tree.query_nearest(geometry), expected)
  1571. @pytest.mark.parametrize(
  1572. "geometry,expected",
  1573. [
  1574. (Point(0.5, 0.5), [0]),
  1575. ([Point(0.5, 0.5)], [[0], [0]]),
  1576. # at junction between 2 lines, will return both
  1577. (Point(2, 2), [1, 2]),
  1578. ([Point(2, 2)], [[0, 0], [1, 2]]),
  1579. # contains one line, intersects with another
  1580. (box(0, 0, 1, 1), [0, 1]),
  1581. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  1582. # overlaps 2 lines
  1583. (box(0.5, 0.5, 1.5, 1.5), [0, 1]),
  1584. ([box(0.5, 0.5, 1.5, 1.5)], [[0, 0], [0, 1]]),
  1585. # second box overlaps 2 lines and intersects endpoints of 2 more
  1586. ([box(0, 0, 0.5, 0.5), box(3, 3, 5, 5)], [[0, 1, 1, 1, 1], [0, 2, 3, 4, 5]]),
  1587. (shapely.buffer(Point(2.5, 2.5), 1), [1, 2, 3]),
  1588. ([shapely.buffer(Point(2.5, 2.5), 1)], [[0, 0, 0], [1, 2, 3]]),
  1589. (shapely.buffer(Point(3, 3), 0.5), [2, 3]),
  1590. ([shapely.buffer(Point(3, 3), 0.5)], [[0, 0], [2, 3]]),
  1591. # multipoints at endpoints of 2 lines each
  1592. (MultiPoint([[5, 5], [7, 7]]), [4, 5, 6, 7]),
  1593. ([MultiPoint([[5, 5], [7, 7]])], [[0, 0, 0, 0], [4, 5, 6, 7]]),
  1594. # second point in multipoint at endpoints of 2 lines
  1595. (MultiPoint([[5.5, 5], [7, 7]]), [6, 7]),
  1596. ([MultiPoint([[5.5, 5], [7, 7]])], [[0, 0], [6, 7]]),
  1597. # multipoints are equidistant from 2 lines
  1598. (MultiPoint([[5, 7], [7, 5]]), [5, 6]),
  1599. ([MultiPoint([[5, 7], [7, 5]])], [[0, 0], [5, 6]]),
  1600. ],
  1601. )
  1602. def test_query_nearest_lines(line_tree, geometry, expected):
  1603. assert_array_equal(line_tree.query_nearest(geometry), expected)
  1604. @pytest.mark.parametrize(
  1605. "geometry,expected",
  1606. [
  1607. (Point(0, 0), [0]),
  1608. ([Point(0, 0)], [[0], [0]]),
  1609. (Point(2, 2), [2]),
  1610. ([Point(2, 2)], [[0], [2]]),
  1611. # 2 polygons in tree overlap point
  1612. (Point(0.5, 0.5), [0, 1]),
  1613. ([Point(0.5, 0.5)], [[0, 0], [0, 1]]),
  1614. # box overlaps multiple polygons
  1615. (box(0, 0, 1, 1), [0, 1]),
  1616. ([box(0, 0, 1, 1)], [[0, 0], [0, 1]]),
  1617. (box(0.5, 0.5, 1.5, 1.5), [0, 1, 2]),
  1618. ([box(0.5, 0.5, 1.5, 1.5)], [[0, 0, 0], [0, 1, 2]]),
  1619. ([box(0, 0, 1, 1), box(3, 3, 5, 5)], [[0, 0, 1, 1, 1], [0, 1, 3, 4, 5]]),
  1620. (shapely.buffer(Point(2.5, 2.5), HALF_UNIT_DIAG), [2, 3]),
  1621. ([shapely.buffer(Point(2.5, 2.5), HALF_UNIT_DIAG)], [[0, 0], [2, 3]]),
  1622. # completely overlaps one polygon, touches 2 others
  1623. (shapely.buffer(Point(3, 3), HALF_UNIT_DIAG), [2, 3, 4]),
  1624. ([shapely.buffer(Point(3, 3), HALF_UNIT_DIAG)], [[0, 0, 0], [2, 3, 4]]),
  1625. # each point in multi point intersects a polygon in tree
  1626. (MultiPoint([[5, 5], [7, 7]]), [5, 7]),
  1627. ([MultiPoint([[5, 5], [7, 7]])], [[0, 0], [5, 7]]),
  1628. (MultiPoint([[5.5, 5], [7, 7]]), [5, 7]),
  1629. ([MultiPoint([[5.5, 5], [7, 7]])], [[0, 0], [5, 7]]),
  1630. (MultiPoint([[5, 7], [7, 5]]), [6]),
  1631. ([MultiPoint([[5, 7], [7, 5]])], [[0], [6]]),
  1632. ],
  1633. )
  1634. def test_query_nearest_polygons(poly_tree, geometry, expected):
  1635. assert_array_equal(poly_tree.query_nearest(geometry), expected)
  1636. @pytest.mark.parametrize(
  1637. "geometry,max_distance,expected",
  1638. [
  1639. # using unset max_distance should return all nearest
  1640. (Point(0.5, 0.5), None, [0, 1]),
  1641. ([Point(0.5, 0.5)], None, [[0, 0], [0, 1]]),
  1642. # using large max_distance should return all nearest
  1643. (Point(0.5, 0.5), 10, [0, 1]),
  1644. ([Point(0.5, 0.5)], 10, [[0, 0], [0, 1]]),
  1645. # using small max_distance should return no results
  1646. (Point(0.5, 0.5), 0.1, []),
  1647. ([Point(0.5, 0.5)], 0.1, [[], []]),
  1648. # using small max_distance should only return results in that distance
  1649. ([Point(0.5, 0.5), Point(0, 0)], 0.1, [[1], [0]]),
  1650. ],
  1651. )
  1652. def test_query_nearest_max_distance(tree, geometry, max_distance, expected):
  1653. assert_array_equal(
  1654. tree.query_nearest(geometry, max_distance=max_distance), expected
  1655. )
  1656. @pytest.mark.parametrize(
  1657. "geometry,max_distance",
  1658. [
  1659. (Point(0.5, 0.5), 0),
  1660. ([Point(0.5, 0.5)], 0),
  1661. (Point(0.5, 0.5), -1),
  1662. ([Point(0.5, 0.5)], -1),
  1663. ],
  1664. )
  1665. def test_query_nearest_invalid_max_distance(tree, geometry, max_distance):
  1666. with pytest.raises(ValueError, match="max_distance must be greater than 0"):
  1667. tree.query_nearest(geometry, max_distance=max_distance)
  1668. def test_query_nearest_nonscalar_max_distance(tree):
  1669. with pytest.raises(ValueError, match="parameter only accepts scalar values"):
  1670. tree.query_nearest(Point(0.5, 0.5), max_distance=[1])
  1671. @pytest.mark.parametrize(
  1672. "geometry,expected",
  1673. [
  1674. (Point(0, 0), ([0], [0.0])),
  1675. ([Point(0, 0)], ([[0], [0]], [0.0])),
  1676. (Point(0.5, 0.5), ([0, 1], [0.7071, 0.7071])),
  1677. ([Point(0.5, 0.5)], ([[0, 0], [0, 1]], [0.7071, 0.7071])),
  1678. (box(0, 0, 1, 1), ([0, 1], [0.0, 0.0])),
  1679. ([box(0, 0, 1, 1)], ([[0, 0], [0, 1]], [0.0, 0.0])),
  1680. ],
  1681. )
  1682. def test_query_nearest_return_distance(tree, geometry, expected):
  1683. expected_indices, expected_dist = expected
  1684. actual_indices, actual_dist = tree.query_nearest(geometry, return_distance=True)
  1685. assert_array_equal(actual_indices, expected_indices)
  1686. assert_array_equal(np.round(actual_dist, 4), expected_dist)
  1687. @pytest.mark.parametrize(
  1688. "geometry,exclusive,expected",
  1689. [
  1690. (Point(1, 1), False, [1]),
  1691. ([Point(1, 1)], False, [[0], [1]]),
  1692. (Point(1, 1), True, [0, 2]),
  1693. ([Point(1, 1)], True, [[0, 0], [0, 2]]),
  1694. ([Point(1, 1), Point(2, 2)], True, [[0, 0, 1, 1], [0, 2, 1, 3]]),
  1695. ],
  1696. )
  1697. def test_query_nearest_exclusive(tree, geometry, exclusive, expected):
  1698. assert_array_equal(tree.query_nearest(geometry, exclusive=exclusive), expected)
  1699. @pytest.mark.parametrize(
  1700. "geometry,expected",
  1701. [
  1702. (Point(1, 1), []),
  1703. ([Point(1, 1)], [[], []]),
  1704. ],
  1705. )
  1706. def test_query_nearest_exclusive_no_results(tree, geometry, expected):
  1707. tree = STRtree([Point(1, 1)])
  1708. assert_array_equal(tree.query_nearest(geometry, exclusive=True), expected)
  1709. @pytest.mark.parametrize(
  1710. "geometry,exclusive",
  1711. [
  1712. (Point(1, 1), "invalid"),
  1713. # non-scalar exclusive parameter not allowed
  1714. (Point(1, 1), ["also invalid"]),
  1715. ([Point(1, 1)], []),
  1716. ([Point(1, 1)], [False]),
  1717. ],
  1718. )
  1719. def test_query_nearest_invalid_exclusive(tree, geometry, exclusive):
  1720. with pytest.raises(ValueError):
  1721. tree.query_nearest(geometry, exclusive=exclusive)
  1722. @pytest.mark.parametrize(
  1723. "geometry,all_matches",
  1724. [
  1725. (Point(1, 1), "invalid"),
  1726. # non-scalar all_matches parameter not allowed
  1727. (Point(1, 1), ["also invalid"]),
  1728. ([Point(1, 1)], []),
  1729. ([Point(1, 1)], [False]),
  1730. ],
  1731. )
  1732. def test_query_nearest_invalid_all_matches(tree, geometry, all_matches):
  1733. with pytest.raises(ValueError):
  1734. tree.query_nearest(geometry, all_matches=all_matches)
  1735. def test_query_nearest_all_matches(tree):
  1736. point = Point(0.5, 0.5)
  1737. assert_array_equal(tree.query_nearest(point, all_matches=True), [0, 1])
  1738. indices = tree.query_nearest(point, all_matches=False)
  1739. # result is dependent on tree traversal order; may vary across test runs
  1740. assert np.array_equal(indices, [0]) or np.array_equal(indices, [1])
  1741. def test_strtree_threaded_query():
  1742. ## Create data
  1743. polygons = shapely.polygons(np.random.randn(1000, 3, 2))
  1744. # needs to be big enough to trigger the segfault
  1745. N = 100_000
  1746. points = shapely.points(4 * np.random.random(N) - 2, 4 * np.random.random(N) - 2)
  1747. ## Slice parts of the arrays -> 4x4 => 16 combinations
  1748. n = int(len(polygons) / 4)
  1749. polygons_parts = [
  1750. polygons[:n],
  1751. polygons[n : 2 * n],
  1752. polygons[2 * n : 3 * n],
  1753. polygons[3 * n :],
  1754. ]
  1755. n = int(len(points) / 4)
  1756. points_parts = [
  1757. points[:n],
  1758. points[n : 2 * n],
  1759. points[2 * n : 3 * n],
  1760. points[3 * n :],
  1761. ]
  1762. ## Creating the trees in advance
  1763. trees = []
  1764. for i in range(4):
  1765. left = points_parts[i]
  1766. tree = STRtree(left)
  1767. trees.append(tree)
  1768. ## The function querying the trees in parallel
  1769. def thread_func(idxs):
  1770. i, j = idxs
  1771. tree = trees[i]
  1772. right = polygons_parts[j]
  1773. return tree.query(right, predicate="contains")
  1774. with ThreadPoolExecutor() as pool:
  1775. list(pool.map(thread_func, itertools.product(range(4), range(4))))