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.

I have a list of points and a list of polygons both in a separate shape files. I am trying to write a ArcPy script to loop through the points, and then the polygons respectively and tell me when a point is within a polygon. I have tried both the contains and touches methods on polygons but haven't had any luck. Any help would be appreciated. So far I have the following:

    pointFeatureClass = "C:\\Data\\POINTS2.dbf"
    pointFieldList = arcpy.ListFields(pointFeatureClass)
    for field in pointFieldList:
        curField = field.name 
        arcpy.AddMessage(curField)

    ptRows = arcpy.SearchCursor("C:\\Data\\POINTS2.dbf", "", "", "CID;Shape","CID;Shape") 
    for ptRow in ptRows: 
        #Now loop through poly's testing for intersect
        shpRows = arcpy.SearchCursor("C:\\Data\\counties.shp", "", "", "NAME;Shape","NAME;Shape")
        for shpRow in shpRows: 
            #arcpy.AddMessage("Checking for point "+ str(ptRow.CID) +" in: " + shpRow.NAME)
            if shpRow.Shape.touches(ptRow.Shape):
                arcpy.AddMessage("Got it: " + shpRow.NAME)
share|improve this question
    
Is there a reason you are coding your own Python method for this, rather than just using Spatial Join? Also, your pointFeatureClass refers to a .dbf file, not a shapefile, so Arcpy will not read the Shape field as geometry. –  dmahr Feb 26 '13 at 16:02
    
My mistake I am also trying this on a dbf that has point x/y coordiantes and converting them to point with arcpy.Point(ptRow.xpos, ptRow.ypos). But, the reason I am trying to script this is that there are 17 million points and I thought that it might be less likely to hang in ArcGIS. –  geo derek Feb 26 '13 at 16:07
2  
Is your polygon in the same coordinate system as the point X/Y coordinates? I just tested it with the contains method, and it worked fine. –  dmahr Feb 26 '13 at 16:25
    
You sir are a gentleman and a scholar. Both of your replies steered me in the correct direction. Everything is working now. Thanks for your help. –  geo derek Feb 26 '13 at 16:46
    
Glad I could help. I've expanded the comments into an answer below, so you can close out the question. –  dmahr Feb 26 '13 at 17:14
add comment

1 Answer 1

up vote 4 down vote accepted

When using ArcPy geometry methods like contains, the coordinates are taken literally with no reprojection. If your XY coordinates are in UTM and your polygon coordinates are in Albers, ArcPy will not correct for the difference, and as a result the contains method will return False.

From a performance perspective, your current "brute force" approach is fairly inefficient: it is looping through all 17 million points times the number of polygons, and performing a computationally "expensive" geometry operation. The built-in Spatial Join tool probably has some space-partitioning algorithms to improve performance, but I think you are right to worry that a Spatial Join with so many points might hang. There are still some other ways you could improve performance:

  • Add a break statement when the contains method returns true. This will make it so your script doesn't continue iterating through other polygons after it has found one that intersects the point.

  • Add some basic partitioning to your points, where possible. Find the extent of your entire polygon layer (not just individual polygon features) and store the XMin, XMax, YMin, and YMax values. If the point has an X coordinate less than the XMin or greater than the XMax, or a Y coordinate less than the YMin or greater than the YMax, then the point must be outside of the polygon layer's extent and can be skipped before even doing any iterating using a continue statement.

  • Apply the same kind of partitioning described in the above bullet to each shape. The contains method of the ArcPy polygon performs a classic "point-in-polygon" search which involves finding the number of intersections between a line going through the point and the perimeter of the shape. This doesn't use any really advanced trig or calc, but it still takes awhile to run, especially for polygons with lots of vertices. You can avoid this by getting the XMin, XMax, YMin, and YMax values of the Extent property of the Polygon. Just as before, if the point has an X coordinate less than the XMin or greater than the XMax, or a Y coordinate less than the YMin or greater than the YMax, then the point must be outside of the polygon extent and can be skipped using a continue statement.

  • If you are feeling brave, you can leave Arc entirely and try using Shapely or another library with less overhead to do your Point-In-Polygon checking. You will want to avoid oddities and edge cases (like multipart polygon features) but I've found that libraries like Shapely have better performance than ArcPy.

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.