[GIS] Improving performance of Clip using OGR

clipgdalogrperformancepython

I have a set of points which I am trying to clip to a set of polygons. The files have many features (150,000 polys, 8,000 points). I'd like to automate this process carry it out once a day, hence I am using python/gdal/ogr to script the operation. However, it is prohibitively slow. Performing the same operation using the Clip tool in ArcGIS gives results in mere minutes! Here is one ogr call I've tried: 

ogr2ogr -clipsrc polygon.shp clipped_output.shp points.shp

and it has been running for 5 hours. I've also tried using a VRT file in combination with an SQL statement to count the number of points falling inside polygons using: 

'ogrinfo -dialect SQLITE -sql "SELECT COUNT(*) FROM polygons, points WHERE ST_Intersects(polygon.geometry,points.geometry)" input.vrt'

This has also taken ~ 5 hours so far on a virtual machine with plenty of RAM and processing power. Are there other techniques or methods that I do not know about that could reduce computation time? Or am I just missing something?

Best Answer

i'm not sure about the ogr2ogr command line tools, though i tend to do similar processes using ogr in python - something like: (assuming these are shapefiles)

from osgeo import ogr

def ogrClip(polys, points, clipppedPoints): 
    # not certain if a separate driver is required for each shapefile?   
    shpdrv = ogr.GetDriverByName('ESRI Shapefile')
    polyDS = shpdrv.Open(polys)
    ptsDS = shpdrv.Open(points)
    polys = polyDS.GetLayer(0)
    pts = ptsD.GetLayer(0)

    clipDS = shpdrv.CreateLayer(clippedPoints, geom_type=ogr.wkbPoint)

    #create desired fields in clipDS - one for example
    fieldDef=ogr.FieldDefn('myIntegerColumn',ogr.OFTInteger)
    clipDS.CreateField(fieldDef)

    # loop through either polys or points and set spatial filter
    # sometimes i reverse this order to test for speed..
    for poly in polys:
        polygeom = poly.GetGeometryRef()
        pts.SetSpatialFilter(polygeom)
        for pt in pts:
            # you could perform intersection here to be certain
            # the points lie within the polygon, but simple geometries
            # seem to work ok without it
            newfeat = ogr.Feature(clipDS.GetLayerDefn())
            newfeat.SetGeometry(pt.GetGeometryRef()) # i think this would work
            newfeat.SetField('myIntegerColumm', pt.GetField('originalColumn'))
            # ....etc.
            clipDS.CreateFeature(newfeat)
            newfeat = None

    shpdrv = None

i didn't test this at the moment - but hopefully you get the general idea.

Related Question