3

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.

3

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

enter image description here

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   

enter image description here

  • My source file is a GDB file and I use it twice to compare two layers... one layer is a surface layer (polygons) and the other layer is points. Should I do anything differently on the load shape files step? – George Jun 7 '16 at 19:46
  • I forgot to mention my source file is a GDB file and I use it twice to compare two layers... one layer is a surface layer (polygons) and the other layer is points. Should I do anything differently on the load shape files step? Other than that I'm using the page you cited as the example... my code is similar to yours up to your step 2... I have step 3's for-loop inside of step 2's for-loop. My "step 3" is more like the page you cited... yours is different, but I will change accordingly to see if I can get past my current impediment. – George Jun 7 '16 at 19:52
  • The example here is with two polygons layers, for a point layer and a polygon layer look at the references cited. – gene Jun 7 '16 at 19:59
  • In that case, my code breaks when trying to pass a 3-D coordinate point to R-Tree... also just wondering is your "if geomb.intersects(geomb):" is a typo... seems like one of the geomb's should be a geomg. – George Jun 7 '16 at 20:00
  • THANK YOU!!! With your recent answer I was able to piece together things to parse through my test GDB in less than 1 minute via a virtual box in spited of my nuanced differences! And that would be going through all water surfaces / polygons. I still need to figure out how to filter for a specific feature (perennial waters) to cut things down even further but I'm not complaining. :) Many thanks again! – George Jun 7 '16 at 21:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.