Vjylku

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.

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]