Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
5  
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
1  
I think you mean ArcInfo workstation. That old beast coupled with AML was impressive! –  OptimizePrime Nov 22 '11 at 16:31
add comment

4 Answers

up vote 3 down vote accepted

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

share|improve this answer
    
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
add comment

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.

share|improve this answer
    
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
add comment

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.

share|improve this answer
    
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
add comment

This is an old post, so this was not available at the time but the Data Access module (v10.1) is super fast and using that and geometry objects allows you to do that type of query faster than the standard toolbox. It is based on NUMPY arrays which is why it is so powerful and fast.

http://resources.arcgis.com/en/help/main/10.1/index.html#/SearchCursor/018w00000011000000/ http://resources.arcgis.com/en/help/main/10.1/index.html#//018z00000070000000

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.