I'm trying to characterize two different python scripts and how fast they rip through the same file GDB on a thin client remotely accessing a Linux server where the FGDB and scripts are stored. The goal is to identify any building found to be contained within perennial waters.

The server I remote into has 16GB of system RAM and an 8GB swap file. The CPU is an 8-core Xeon E5-2695 v2 @ 2.40GHz.

In one approach, I use Shapely, R-Tree, and Fiona to reach the same results as PostgreSQL v9.3.14, OGR, Psycopg, and PostGIS. The same server hosts both approaches.

Long story short, my PostGIS approach spanks my Shapely approach. I feel there's room to close the gap with my Shapely approach to be more "competitive" to my PostGIS approach.

I'm using 3 different FGDBs.

  • FGDB 1 is about 82.7MB in size with 6,070,025 "raw" permutations
  • FGDB 2 is about 612.8MB in size with 205,425,792 "raw" permutations
  • FGDB 3 is about 11.4GB in size with 104,055,912,442 "raw" permutations

Shapely's approach yields (read into mem time, analysis time, total time):

  • FGDB 1: 19.71 sec, 141.71 sec, 161.42 sec
  • FGDB 2: 272.56 sec, 145.3 sec, 417.86 sec
  • FGDB 3: 1,313.73 sec, 4,493 sec, 5,806.73 sec

PostGIS's approach yields (conversion time, analysis time, total time):

  • FGDB 1: 12 sec, 8.13 sec, 20.13 sec
  • FGDB 2: 58.04 sec, 28.14 sec, 86.18 sec
  • FGDB 3: 1,268.14 sec, 397.43 sec, 1,665.57 sec

PostGIS is 8x, 4.8x, and 3.5x faster than Shapely respectively for each FGDB.

My Shapely code is below:

for bldg in buildingLayerList:
    with fiona.open(sourceDatabaseFile, 'r', layer=bldg) as sourceFile1:
        for building in sourceFile1:
            buildingFeatures.append((building['id'], building['geometry']))

    #Use current list of buildings and compare against all hydrography layers
    for hydro in hydrographyLayerList:
        with fiona.open(sourceDatabaseFile, 'r', layer=hydro) as sourceFile2:
            for water in sourceFile2:
                waterFeatures.append((water['id'], water['geometry'], water['properties']))

        #from rtree import index
        idx = index.Index() #Create an R-Tree index and store the features in it (bounding box)
        for pos, poly in enumerate(waterFeatures):
            idx.insert(pos, shape(poly[1]).bounds)

        for building in buildingFeatures:
            buildingShape = shape(building[1])

            #if the geometry intersects the spatial index
            for id in idx.intersection(buildingShape.bounds):
                water = waterFeatures[id]

                #Check to see if current hydrography feature intersection is a perennial water
                if (('ZI024_HYP' in water[2]) and (water[2]['ZI024_HYP'] == 1)):
                    waterShapeGeometry = shape(water[1])
                    numBoundedBldgHydroFeaturesEvaluated += 1

                    #Check to see if perennial water actually contains building feature
                    if waterShapeGeometry.contains(buildingShape):

            if len(discoveredBuildingsContainedInWaters) > 0:
                totalBldgIDs.update({bldg + " vs " + hydro: sorted(discoveredBuildingsContainedInWaters)})

        discoveredBuildingsContainedInWaters[:] = [] # Just in case
        waterFeatures[:] = [] # Just in case
    buildingFeatures[:] = [] # Just in case

Am I incorrectly implementing a spatial index / bounding box or am I inefficiently reading in every feature into memory?


I rewrote the block of my code to avoid unnecessarily visiting the water layer twice.

idx = index.Index() #Create an R-Tree index and store the features in it (bounding box)
for pos, poly in enumerate(waterFeatures):
    idx.insert(pos, shape(poly[1]).bounds)

for i,bldgs in enumerate(buildingFeatures):
    building = shape(bldgs[1])

    pt2D = (list(building.coords)[0][:2])
    for j in idx.intersection(pt2D):
        if building.within(shape(waterFeatures[j][1])):

The downside is these changes have little to change performance speed.

  • 2
    Unless you have some significantly more complex files, I'd suggest that any optimisation of either approach is a waste of time. – BradHards Nov 17 '16 at 21:05
  • "GDB" seems to be used as a generic term to reference spatial data. This is not a conventional use of that term, which adds confusion the question. – Vince Nov 18 '16 at 0:18
  • 1
    You should profile your code to see which parts are taking the most time and focus on those. The two things that jump out to me are that you're calling shape() on the water features twice, and that you're not using the bulk import to rtree using a generator. – Snorfalorpagus Nov 18 '16 at 11:20
  • @Vince When I stated GDB, I'm meant to say ESRI-generated GDB "file" folders. Should I have prefaced that in my original post? – George Nov 18 '16 at 17:14
  • @Snorfalorpagus, I'm learning Shapely as I go. With regard to calling Shapely for the water features twice, how would you recommend I reduce that down to one call? My first call is to convert a water geometry into a Shapely shape. The next call is to do an actual evaluation. Are you suggesting that Shapely can do the latter without the former? – George Nov 18 '16 at 17:22

Your Answer

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

Browse other questions tagged or ask your own question.