[GIS] Efficient arcpy method of points in polygon

arcgis-10.0arcpypython

Coming from a background of Mapinfo and Manifold, I frequently find myself frustrated with the ArcGIS SQL toolset. In both these former programs it is trivial to write a query asking, for instance, how many points fall within each polygon. This can be filtered using where clauses, e.g. all the points of type a within each polygon and then all tbe points of type b within each polygon. Since the query result is an in-memory table, it is quick and uses up very little disk space.

By contrast, Spatial Join in ArcGIS seems quite slow, possibly because of disk IO, uses up more disk space than necessary, and I personally find it cumbersome.

I know there are alternatives, Hawths Tools for instance, but I was trying to do something in arcpy. My first thought was to loop through each geometry in both the polygon set and point set and test whether the point was contained by the polygon. This would then enable all sorts of things to be calculated. I wrote a test that just looped through all the points for one polygon:

import arcpy
arcpy.env.workspace = r"path\to\workspace"
print "Start"
geom = arcpy.Point(530000,180000)
cursor = arcpy.SearchCursor("FCName")
for rowid in cursor:
     poly = rowid.Shape
     if poly.contains(geom):
          print rowid.OBJECTID

However, this takes a few seconds for each polygon, so with a dataset size of 4,000 polygons, I don't think this is going to work particularly quickly.

Does anyone have any ideas for making this more efficient, or a function that I'm missing?

Best Answer

Would using the Intersect tool work for you? Each point in the resulting feature class would also have the attributes of the polygon it is contained in, and you could write a script that loops through those points and calculates whatever statistics you need for each polygon ID (min, max, mean, sum, count, etc.).

Related Question