1

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?

  • 1
    Does your first command improve if you use the SQL statement in the -clipsrcsql flag? Also try ST_Within instead of ST_Intersects. It also might help adding the -dialect SQLITE option to your first command (like in your second command), if possible. – SaultDon Nov 12 '15 at 23:17
  • 1
    The first query might be faster if you make an union of the 150000 polygons to be used as -clipsrs. For the second problem: do not count how many points are inside each of the 150000 polygons. Make a query that resolves within which polygon each point is. That makes only 8000 computations, and for optimizing those I suggest to read this document gaia-gis.it/spatialite-3.0.0-BETA1/WorldBorders.pdf. – user30184 Nov 13 '15 at 0:31
  • @SaultDon when I use the -clipsrcsql flag the command returns the same input points almost instantly, but they are not clipped to the polygons. Adding -dialect SQLITE gives the same result. Adding -dialect SQLITE and using -clipsrc (as opposed to -clipsrcsql) may do the trick - I'll see in a few hours but the filesize on clipped_output.shp does not appear to be changing which makes me think it will not work. I'm also testing out the ST_Within option on the second command right now. – geoeye Nov 13 '15 at 15:42
  • @user30184 thanks for your suggestions. I looked at the link you provided but I do not think it's very useful in my case as my polygons are very simple geometries, most containing just 4 lat/lon coordinates. As for a query that resolves which polygon each point is in, here's my attempt: ogrinfo -dialect SQLITE -sql "SELECT COUNT(*) FROM points, polygons WHERE ST_WITHIN(points.geometry, polygons.geometry)" input.vrt I incorporated SaultDon 's suggestion to use ST_Within instead of ST_Intersects. Is that method likely to work, and does my query look to be reasonable? It is currently running – geoeye Nov 13 '15 at 15:54
  • @User847387435, did it work??? I need the same! – J Kelly Jun 1 '16 at 15:23
1

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.

0

If your clipping zone is small considering the whole vector area, you could redice a lot the time using the --spat xmin ymin xmax ymax. The only problem is that you need to calculate xmin ymin xmax ymax before.

It will be something like

ogr2ogr --spat xmin ymin xmax ymax -clipsrc polygon.shp clipped_output.shp points.shp

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.