How To Properly Implement a Bounding Box for Shapely & Fiona?

Question

I am trying to spatially analyze two layers within same GDB file using Python / Shapely.

Another user was kind enough to give me a procedure to go by in trying to reduce the number of senseless Shapley "within" calls by doing the following:

1. Read the features in your 1st layer into a list like list(gdb)
2. For each feat in this list do a list(gdbStructures.items(bbox=fiona.bounds(feat))) to get a list of the features in the second dataset that have bounding boxes intersecting the bounding box of feat.
3. For each element of that result list, make a a.within(b) test as you were doing.

Shapely's within() can be expensive if your shapes have a lot of detail. The approach above will use the spatial indexes in your GDB files to make those calculations only when needed.

Below is my code snippet:

with fiona.drivers():
    try:
        if os.path.isdir(sourceDatabaseFile):

            numBLDGwithinPerennialWaters = 0
            listOfPerennialWaters = []
            listOfAllBuildings = []

            with fiona.open(sourceDatabaseFile, 'r', layer="HydrographySrf") as gdb: # Pointer to HydrographySrf layer of source file
                with fiona.open(sourceDatabaseFile, 'r', layer="StructurePnt") as gdbStructures: # Pointer to StructurePnt layer of source file

                    for hydroFeature in gdb:
                        if (('ZI024_HYP' in hydroFeature['properties']) and (hydroFeature['properties']['ZI024_HYP'] > 0)):
                                listOfPerennialWaters.append(hydroFeature)

                    for waterFeature in listOfPerennialWaters:
                        for building in gdbStructures:
                            list(gdbStructures.items(bbox=fiona.bounds(waterFeature))) 
                            if shape(waterFeature['geometry']).contains(shape(building['geometry'])):

                                ...

I have tried to do variations of the list command where the bbox is invoked... such as assigning ListofAllBuildings to whatever is returned from the list command or appending ListOfAllBuildings. No matter what I seem to do I get zero after the looping finishes (presumably storing on the last List command called).

I do know at some iterations of the loop where I call the list command, I do get non-zero lengths.

Answer 1

I show you here a simple example of the use of bounding boxes in a spatial index (intersection of two polygon shapefiles)

1) load the shapefiles

import fiona
from shapely.geometry import shape
blue = fiona.open("poly1.shp")
green = fiona.open("poly2.shp")

2) create a spatial index with the bounding boxes of the blue layer (you can use rtree, pyrtree (R-tree indexes) or Pyqtree (Quadtree spatial index) or ....)

from rtree import index
idx = index.Index()
for pos, poly in enumerate(blue):
   idx.insert(pos, shape(poly['geometry']).bounds)

3) Now with the bounding spatial index, you iterate only through the geometries which have a chance to intersect with your current geometry.

for featg in green:
   geomg = shape(featg['geometry'])
   # if the geometry intersects the spatial index
   for id in idx.intersection(geomg.bounds):
          featb = blue[id]
          geomb= shape(featb['geometry'])
          if geomb.intersects(geomb):
               print '{} intersects {}'.format(geomg['id'], geomb['id'])
  0 intersects 0
  0 intersects 2
  1 intersects 3

And you can compute the intersections.

A faster solution was proposed by Mike T in Faster way of polygon intersection with shapely

from shapely.ops import unary_union
for poly in green:
    merged_cells = unary_union([shape(blue[pos]['geometry']) for pos in idx.intersection(shape(poly['geometry']).bounds)])
    shape(poly['geometry']).intersection(merged_cells).wkt