Hfrhyu

I have a set of Polygons (in Shapely) that I wish to find all intersections that are comprised of overlaps of two or more polygons. I have been working based off the code found in this answer. I have a code so far that works in the if 1: statement below, but it takes a very long time to run. The code that is in the if 0: statement, runs quickly, but produces an output with some errors. See below for an example of the error.

from shapely.geometry import Point, LineString, shape, mapping
from shapely.ops import cascaded_union
import fiona
import shapely
from shapely.ops import unary_union, polygonize
from osgeo import ogr
from shapely.strtree import STRtree

with fiona.open("Untitled.shp") as layer:
 subset=layer
 spoof=[]
 for n,pol in enumerate(subset):
 if type(shape(pol['geometry'])) == shapely.geometry.polygon.Polygon:
 spoof.append(shape(pol['geometry']))
 elif type(shape(pol['geometry'])) == shapely.geometry.multipolygon.MultiPolygon:
 print('theres one here')
 tmp=list(shape(pol['geometry']))
 for x in tmp:
 spoof.append(x)

 rings = [LineString(list(pol.exterior.coords)) for pol in spoof]
 union = unary_union(rings)

 tree = STRtree(spoof)
 result = [geom for geom in polygonize(union)]

 final=[]
 if 1:
 for k,item1 in enumerate(result):
 count=0
 for n,item2 in enumerate(spoof):
 if item1.intersection(item2).area >1e-8:
 count=count+1
 if count>=2:
 final.append(item1)

 if 0:
 for k,item in enumerate(result):
 touch=tree.query(item.centroid)
 if len(touch)>=2:
 final.append(item)

Below, the orange represents all the shapes that I would like to check if they are produced (by unary union in Shapely) by two or more overlaps, the cyan represents the output of the if-statement that works (if 1 in above code), and the pink represents the incorrect selected polygons by the if-statement that doesn't work (if 0 in above code). The top left corners of the text boxes in the pictures represent that point being selected, and as you can see the second picture, the pink area only has one polygon from the STEREO_OBSERVATIONS (which is the set of polygons that are being checked for overlap) so it should not be selected. By contrast, the left adjacent cyan square to the pink square highlighted in the third picture has two polygons from STEREO_OBSERVATIONS and as such is correctly selected by the cyan result. The layering of the polygon layers is shown in the first picture.

What is the best way to query the STRtree structure so as to check to intersection or overlap? I have tried to use the centroid in the code, but that hasn't seemed to produce the correct output.

Fig. 1 (Legend/layering)

Fig 2. (Incorrect selection by fast algorithm [if 0 in code] in pink)

Fig 3. (Correct selection by slow algorithm [if 1 in code]

The solution ended up being to change the if-statement to:

 if 1:
 for k,item in enumerate(result):
 print(len(result),k)
 touch=tree.query(item.centroid)
 inside=[]
 for x in touch:
 if x.contains(item.centroid):
 inside.append(item)
 if len(inside)>=2:
 final.append(item)

The reason for the original error was that the tree.query(<>) only returns the polygons that are likely to interact with the query <> object based on the bounds (top, bottom, left, right) of the objects inside of the tree. The solution was to check whether or not the centroid of the polygon was contained within each member of the list generated by the query. In some cases, the objects returned by the query were in the vicinity of the query object but did not in fact contain it.