topic arcpy points in polygon check in Python Questions

arcpy points in polygon check

JohnMorgan — Tue, 26 Feb 2013 13:51:39 GMT

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)

Re: arcpy points in polygon check

MathewCoyle — Tue, 26 Feb 2013 14:08:19 GMT

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

Re: arcpy points in polygon check

JohnMorgan — Sun, 12 Dec 2021 08:08:23 GMT

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

Re: arcpy points in polygon check

MathewCoyle — Tue, 26 Feb 2013 14:50:16 GMT

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

Re: arcpy points in polygon check

JohnMorgan — Thu, 28 Feb 2013 00:25:02 GMT

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.

Re: arcpy points in polygon check

NettoChad — Wed, 27 Aug 2014 22:13:18 GMT

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)