arcpy points in polygon check

7273
5
Jump to solution
02-26-2013 05:51 AM
JohnMorgan
Occasional Contributor
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)
Tags (2)
0 Kudos
1 Solution

Accepted Solutions
MathewCoyle
Frequent Contributor
Is there some reason Spatial Join doesn't work for you?

View solution in original post

5 Replies
MathewCoyle
Frequent Contributor
Is there some reason Spatial Join doesn't work for you?
JohnMorgan
Occasional Contributor
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
0 Kudos
MathewCoyle
Frequent Contributor
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
0 Kudos
NettoChad
New Contributor III

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)

0 Kudos
JohnMorgan
Occasional Contributor
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.
0 Kudos