arcpy points in polygon check

[JohnMorgan]

Hello,
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:

import arcpy  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)

[MathewCoyle · Accepted solution]

Is there some reason Spatial Join doesn't work for you?

[JohnMorgan]

Thanks for you response. I am trying that now.  I should also point out that I am starting with a list of points in a dbf and then converting them to point geometry with

point = arcpy.Point(ptRow.xpos, ptRow.ypos)
 ptGeometry = arcpy.PointGeometry(point) 

And that I am dealing with about 17 million points and 3000 polygons.  Therefore, I was assuming that this would be something to script.  But if the join works then that would be great.

Derek

[MathewCoyle]

That is quite a few points to get through but should be doable. Is this something you will be doing with some frequency? If so you would definitely want to write a script for it, however basic. If it is something you want to optimize and/or will be doing quite often you can get right down to the math and bypass arcpy completely (albeit reinventing the wheel a bit). It follows the concept of the Jordan curve theorem to determine whether a point is inside or outside a polygon. You can utilize the code for the even-odd rule here to get started. http://en.wikipedia.org/wiki/Even-odd_rule

[NettoChad]

I came across this method some time ago. This method does not use arcpy and is much faster. The only issue is that you have to import the geometry of your polygons.

GeospatialPython.com: Point in Polygon 2: Walking the line

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "Point " + str(x) + "," + str(y) + " is IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "Point " + str(x) + "," + str(y) + " is IN"

   n = len(poly)
   inside = False

   p1x,p1y = poly[0]
   for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
         if y <= max(p1y,p2y):
            if x <= max(p1x,p2x):
               if p1y != p2y:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "Point " + str(x) + "," + str(y) + " is IN"
   else: return "Point " + str(x) + "," + str(y) + " is OUT"

# Test a vertex for inclusion
poligono = [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat= -33.416032
lon= -70.593016
print point_in_poly(lat, lon, poligono)
# test a boundary point for inclusion
poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)]
x = 10
y = 1
print point_in_poly(x, y, poly2)

[JohnMorgan]

That is quite a few points to get through but should be doable. Is this something you will be doing with some frequency? If so you would definitely want to write a script for it, however basic. If it is something you want to optimize and/or will be doing quite often you can get right down to the math and bypass arcpy completely (albeit reinventing the wheel a bit). It follows the concept of the Jordan curve theorem to determine whether a point is inside or outside a polygon. You can utilize the code for the even-odd rule here to get started. http://en.wikipedia.org/wiki/Even-odd_rule

Thanks Mathew.  That is interesting reading.  The spatial join did end up working for me best though as you first replied.