Correct syntax for using arcpy.Geometry.touches or .intersects methods

[StevenLambert]

Am trying to test new polygon created in script for overlap against an existing polygon feature class that may contain one or more multippart polygons.  Can seem to get the right syntax.  Keep getting failure on second geometry parameter.....

# First create geometry object to receive existing polygon features

inML = arcpy.geometry()

# Copy in features from fc to geometry object

arcpy.CopyFeatures_management(infc, inML)

# Have polygon footprint in an array, aFoot, create geometry object

test_poly = arcpy.Geometry("polygon", aFoot)

# test for overlap

intersect = inML.touches(test_poly)

or

intersect = inML.intersects(test_poly)

Both fail with a second geometry error....

Any ideas?

Thanks.

[DavidWynne]

Hi Steven,
First, the geometry object you have specified as output doesn't get populated with the geometry as your code is expecting.  That geometry object is just a special way of signifying you want geometry as output; the actual output is really on the left-side.  Something like below would work:

gL = arcpy.CopyFeatures_management(infc, arcpy.Geometry())

In this case, the variable gL is a list of geometries.  If you have just one feature you could then use a relational operator like this:

intersect = gL[0].touches(test_poly)

One small aside.  When creating your Polygon geometry, you could simplify you code like this:

test_poly = arcpy.Polygon(aFoot)

instead of:

test_poly = arcpy.Geometry("polygon", aFoot)

Hope this helps,
-Dave

[StevenLambert]

Thanks. David. 

One additional question,  for multipart geometries (e.g. geometry.isMultipart = True) , do I have to manually cycle through the individual parts of each geometry in the gL to test individual elements for intersection or does the geometry.touches  automatically test all parts each geometry in the gL?

[DavidWynne]

One additional question,  for multipart geometries (e.g. geometry.isMultipart = True) , do I have to manually cycle through the individual parts of each geometry in the gL to test individual elements for intersection or does the geometry.touches  automatically test all parts each geometry in the gL?

No, you won't have to cycle through the parts if working with a multipart geometry.  touches will work against the entire feature.

-Dave

[StevenLambert]

That's great.  I appreciate the help.  Thanks again.

[ChrisMathers]

Thats awesome. If you have two two features, each in a separate feature class, aFeature and bFeature, could this method of checking polygon intersection be used to see if aFeature intersects aFeature and if so write an attribute of bFeature to a field in aFeature?

[DavidWynne]

Thats awesome. If you have two two features, each in a separate feature class, aFeature and bFeature, could this method of checking polygon intersection be used to see if aFeature intersects aFeature and if so write an attribute of bFeature to a field in aFeature?

Hi clm42,
Yes, and no.  You could definitely use this general approach, but geometry objects don't hold attribute information.  So instead of using the approach of getting geometry as output from geoprocessing tools (as most of the code in this thread shows), you'd have to open some cursors.  The cursor will then give you both the geometry and attribute information for a feature.

-Dave

[StevenLambert]

Was not getting convincing and consisten results, so I tried tests in code block.  I think that the "crosses" and "touches" are either wrong or the definitions in the documentation are inconsistent.  What do you think?

#
# test of the various geometry or poly interaction methods
#
import arcpy
import sys

# define coordinate lists

clist1 = [[0,0], [0,10],[10,10],[10,0],[0,0]]
clist2 = [[10,0],[10,10],[20,10],[20,0],[10,0]]

# define geometry objects
#
#   (0, 10) -------(10, 10)---------(20, 10)
#           |  P1  |       |   P2  |
#   (0,  0) -------(10,  0)---------(20,  0)
#
# Intersection along the line (10,0) to (10,10)

aPt = arcpy.Point()
anArray = arcpy.Array()

# Make first polygon

for coordpair in clist1:
    aPt.X, aPt.Y = coordpair[0], coordpair[1]
    anArray.add(aPt)

polygon1 = arcpy.Polygon(anArray)

# Clear array and push in second coordinate list

anArray.removeAll()

# Make second polygon

for coordpair in clist2:
    aPt.X, aPt.Y = coordpair[0], coordpair[1]
    anArray.add(aPt)

polygon2 = arcpy.Polygon(anArray)

# Method Explanation
#
#   contains (second_geometry) Indicates if this geometry contains the other geometry. 
#   crosses (second_geometry) Indicates if the two geometries intersect in a geometry of lesser dimension. 
#   disjoint (second_geometry) Indicates if the two geometries share no points in common.  
#   equals (second_geometry) Indicates if the two geometries are of the same type and define the same set of points in the plane. 
#   overlaps (second_geometry) Indicates if the intersection of the two geometries has the same dimension as one of the input geometries. 
#   touches (second_geometry) Indicates if the boundaries of the geometries intersect. 
#   within (second_geometry) Indicates if this geometry is contained within another geometry. 

# Based on coordinate lists and definitions

# Polygon1 contains Polygon2 = False
# Polygon1 crosses Polygon2 = True Two two-dimensional polygons share a one-dimensional boundary, the line from (10,0) to (10,10)
# Polygon1 disjoint Polygon2 = False
# Polygon1 equals Polygon2 = False
# Polygon1 overlaps Polygon2 = False (Intersection is one-dimensional line (10,0) to (10,10)
# Polygon1 touches Polygon2 = True (Intersect along line (10,0) to (10,10)
# Polygon1 within Polygon2 = False

# Evaluate various polygon methods

test1 = polygon1.contains(polygon2)
test2 = polygon1.crosses(polygon2)
test3 = polygon1.disjoint(polygon2)
test4 = polygon1.equals(polygon2)
test5 = polygon1.overlaps(polygon2)
test6 = polygon1.touches(polygon2)
test7 = polygon1.within(polygon2)

# Evaluate obverse

test8 = polygon2.contains(polygon1)
test9 = polygon2.crosses(polygon1)
test10 = polygon2.disjoint(polygon1)
test11 = polygon2.equals(polygon1)
test12 = polygon2.overlaps(polygon1)
test13 = polygon2.touches(polygon1)
test14 = polygon2.within(polygon1)

print 'Test1  = ', test1
print 'Test2  = ', test2
print 'Test3  = ', test3
print 'Test4  = ', test4
print 'Test5  = ', test5
print 'Test6  = ', test6
print 'Test7  = ', test7

print ' '

print 'Test8  = ', test8
print 'Test9  = ', test9
print 'Test10 = ', test10
print 'Test11 = ', test11
print 'Test12 = ', test12
print 'Test13 = ', test12
print 'Test14 = ', test12


print ' '

[StevenLambert]

I retract concern about crosses test.  Was focusing on planar polygons.  Only way two polgons can cross at a lesser dimesionality is to share only a point e.g. polygons (and lines!) are 2D points are 1D...

Still not sure about touches test though...

[DanPatterson_Retired]

I would test by creating 2 known shapes and see if the test passes or fails.  Recently
http://forums.arcgis.com/threads/9763-Errors-in-arcpy-s-Polygon-class
several issues have been raised about the arcpy polygon and polyline class accuracies which affect the extent of features as well as other parameters.  To test create two adjacent shapes see if they intersect, offset them by a finite epsilon amount and see if the test fails.