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)
Solved! Go to Solution.
point = arcpy.Point(ptRow.xpos, ptRow.ypos) ptGeometry = arcpy.PointGeometry(point)
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)
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