# Efficient arcpy method of points in polygon

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?

-
Lightweight? Fast? Simple? ESRI product? It was called ArcView 3.x :-). (Sure, it had its faults, but in many ways ArcGIS makes the old ArcView look really good...) –  whuber Nov 22 '11 at 15:43
I think you mean ArcInfo workstation. That old beast coupled with AML was impressive! –  OptimizePrime Nov 22 '11 at 16:31

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.).

-
this is a good idea, and I think it will be much quicker. Will try it tommorrow –  Stev_k Nov 22 '11 at 17:29
It worked, and was much quicker, thanks –  Stev_k Nov 23 '11 at 11:38
Glad to help :) –  nmpeterson Nov 23 '11 at 15:01

Can you call the spatial join tool from your arcpy script? Maybe write a wrapper method to handle where the output goes and deal with deleting it afterwards. That will probably be quicker than looping through features directly in arcpy.

Otherwise, I'm afraid arcpy is slow. You could use ctypes to call a dll that does the point-in-polygon function for you, but that's a lot more implementation work, and you still have the overhead of reading features in the first place.

-
yes, I could but I don't really like the Spatial Join tool - I want something more lightweight and flexible –  Stev_k Nov 22 '11 at 11:55

I am like you. I am required to use ArcPy and in the current environment I'm working in, using the SpatialJoin or other ArcGIS tools is cumbersome.

Due to other constraints, I can not use this library. Have a look at Shapely which you can download using easy_install I believe.

It has a range of speedy methods to perform the most common and basic of functions. Couple this with shplib and you'll have yourself a fast and light weight alternative.

-
I would love to try Shapely, but unfortunately my constraints (locked down system) means it will be difficult. May be able to have a play though at some point, thanks –  Stev_k Nov 23 '11 at 11:39
Then you may have to look at cooking your own solution. I have a couple of methods that I have in python that do what you're asking. They might be a good start –  OptimizePrime Nov 24 '11 at 20:22