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:
- Read the features in your 1st layer into a list like list(gdb)
- 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.
- 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.