Need help using geometry to test whether point is within a polygon

3129
4
Jump to solution
03-16-2016 03:21 PM
RebeccaStrauch__GISP
MVP Esteemed Contributor

Back in Avenue days, I was able to use something like

if ((checkPoly.Contains(testline).Not)) then  ...

To test wither my line segment (testline) was contained within the study polygon (checkpoly).  I'm now trying to do the sam type of test use the .within method.

I have created my line geometry and polygon geometry, but I'm receiving  Runtime error.  I have also tried it with a pt geometry.  All three geometries are in the same spatial reference and should return a "true" in both cases.  As powerful as geometries are and as often as I use them, I never seem to be able to get the "light bulb" to go off in my head....so any suggestions are welcome.  I have looked at other threads including https://community.esri.com/thread/69044#comment-400410 and may pursue that, but my guess is it is something simple I am missing.

>>> polyGeom = arcpy.CopyFeatures_management(inStudy, arcpy.Geometry())
... ptGeom = arcpy.PointGeometry(arcpy.Point(ptX, ptY, 0, 0, ptid))
... tmpTrans = geo.segmentAlongLine(startDist, endDist)
... 
... print("check for pt within poly")
... ptGeom.within(polyGeom)
... print("check for line within poly")
... tmpTrans.within(polyGeom)
... 
check for pt within poly
Runtime error 
Traceback (most recent call last):
  File "<string>", line 6, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\arcobjects\arcobjects.py", line 771, in within
    return convertArcObjectToPythonObject(self._arc_object.Within(*gp_fixargs([second_geometry])))
ValueError: [<geoprocessing describe geometry object object at 0x1E43D4E0>]
>>> print("check for line within poly")
... tmpTrans.within(polyGeom)
... 
check for line within poly
Runtime error 
Traceback (most recent call last):
  File "<string>", line 2, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\arcobjects\arcobjects.py", line 771, in within
    return convertArcObjectToPythonObject(self._arc_object.Within(*gp_fixargs([second_geometry])))
ValueError: [<geoprocessing describe geometry object object at 0x1E43D4E0>]
>>> 

Sorry about the formatting.....I was trying to include the errors.  Thanks.

tagging python snippets​  in case anyone has any snippets already to do this.

0 Kudos
1 Solution

Accepted Solutions
DarrenWiens2
MVP Honored Contributor

Untested, but it looks like polyGeom is a list of polygons, not a single polygon geometry. Try:

ptGeom.within(polyGeom[0])

edit: now tested with two points, one inside, one outside a single polygon

>>> polyGeom = arcpy.CopyFeatures_management("camps",arcpy.Geometry())
... print polyGeom
... print type(polyGeom)
... with arcpy.da.SearchCursor("new_point",'SHAPE@') as cursor:
...     for row in cursor:
...         print row[0].within(polyGeom[0])
...         
[<Polygon object at 0x6ee1090[0x6ee1380]>]
<type 'list'>
True
False

View solution in original post

4 Replies
DarrenWiens2
MVP Honored Contributor

Untested, but it looks like polyGeom is a list of polygons, not a single polygon geometry. Try:

ptGeom.within(polyGeom[0])

edit: now tested with two points, one inside, one outside a single polygon

>>> polyGeom = arcpy.CopyFeatures_management("camps",arcpy.Geometry())
... print polyGeom
... print type(polyGeom)
... with arcpy.da.SearchCursor("new_point",'SHAPE@') as cursor:
...     for row in cursor:
...         print row[0].within(polyGeom[0])
...         
[<Polygon object at 0x6ee1090[0x6ee1380]>]
<type 'list'>
True
False
RebeccaStrauch__GISP
MVP Esteemed Contributor

hmmm.  that did get rid of the error by my result is false which doesn't make too  much sense.  btw polyGeom did just have one poly in this case.

I may have to look at the input a little bit more.  because the actual spatial reference isn't specified, maybe that is an issue even if the coordinates are just coordinates?

0 Kudos
DarrenWiens2
MVP Honored Contributor

If everything uses the same coordinate system, it doesn't matter (like in my example above), but if they use different CRS or you need to make measurements in some way (e.g. distance along line), then you must specify.

Here are some ways you can get and set CRS for geometry objects:

>>> polys = arcpy.CopyFeatures_management("camps",arcpy.Geometry())
... polyGeom = polys[0]
... poly_sr = polyGeom.spatialReference # get CRS from a geometry
... print poly_sr.name
... points = "proj_point"
... points_sr = arcpy.Describe(points).spatialReference # get CRS from a feature layer
... print points_sr.name
... with arcpy.da.SearchCursor(points,'SHAPE@',spatial_reference=points_sr) as cursor: # set CRS for cursor geometries
...     for row in cursor:
...         print row[0].within(polyGeom.projectAs(points_sr)) # project geometry to match cursor CRS
...         
NAD_1983_UTM_Zone_10N
NAD_1983_BC_Environment_Albers
True
False
RebeccaStrauch__GISP
MVP Esteemed Contributor

Thanks Darren, you example helped.  In this case my in tmpTrans actually has two parts and one part was in, one part just happened to be out...didn't think about that at the time, so I'll will have some additional testing.  In short, given a random point, grab the contour...grab a transect with the pt as the mid-pont, check that it is within the poly...adjust as necessary.  As fate would have it, the random point I selected just happened to be near the top of a hill so part of the contour is a circle...the rest is outside the study.  I had tests for all this in my Avenue program.... 

This is what I ended up with

>>> for pnt in tmpTrans: 
...  print pnt[0].within(polyGeom[0]) 
...  print pnt[0]
... 
False
402044.161016015 1527045.9042345 NaN 1094633.70727898
True
341396.149335118 1561394.85163668 NaN 1097347.1678
>>>

Thanks for your help!

0 Kudos