Get access to polygon containing point using Python arcpy

[Montgomery_CountyPlanning_Depa]

Hello all,

I have a parcels layer in SDE GDB which contains over 300,000 polygons. I want to print attributes of a parcel polygon which "contains" an input point. The following code works but it is slow (takes about a minute to run). Is there a way to re-write this code (possibly using an Intersect or Spatial Join approach) for faster performance? Please provide some sample code. Thanks!

import arcpy

# Input point coordinates
locX = 1303479.87779
locY = 485423.249288


dataMxdPath = r"C:\Projects\Parcels.mxd"
dataMxd = arcpy.mapping.MapDocument(dataMxdPath)
parcelsLayer = arcpy.mapping.ListLayers(dataMxd, "Parcels")[0]

point = arcpy.Point(locX, locY)

for row in arcpy.SearchCursor(parcelsLayer):
    polygonGeom = row.SHAPE
    if polygonGeom.contains(point):
        print row.OBJECTID
        break

[IanMurray · Accepted solution]

Have you tried using a data access cursor(arcpy.da.SearchCursor)?  Its significantly faster than the old cursor.  Also, your current cursor uses all fields, which takes mroe time, you can specify only the geometry field, which should cut down processing time.  I don't see how using any geoprocessing tools would speed things up.

for arcpy.da.SearchCursor usage and syntax:

ArcGIS Help 10.1‌

import arcpy  


# Input point coordinates  
locX = 1303479.87779  
locY = 485423.249288  
  
  
dataMxdPath = r"C:\Projects\Parcels.mxd"  
dataMxd = arcpy.mapping.MapDocument(dataMxdPath)  
parcelsLayer = arcpy.mapping.ListLayers(dataMxd, "Parcels")[0]  


point = arcpy.Point(locX, locY)  


with arcpy.da.SearchCursor(parcelsLayer, ["SHAPE@"] ) as cursor: 


    for row in cursor:
        polygonGeom = row[0]  
        if polygonGeom.contains(point):  
            print row.OBJECTID  
            break  

Also I added the with statement for memory's sake.  Your original code doesn't delete the cursor object after the script ends, but the with takes care of that.

[Montgomery_CountyPlanning_Depa]

@Ian

Thank you so much for your reply! It was helpful. I had to modify your code slightly to get it to print. It works much faster now (takes about 10 seconds instead of taking a minute like before). Modified working code is pasted below:

import arcpy

# Input point coordinates
locX = 1303479.87779
locY = 485423.249288

dataMxdPath = r"C:\Projects\Parcels.mxd"
dataMxd = arcpy.mapping.MapDocument(dataMxdPath)
parcelsLayer = arcpy.mapping.ListLayers(dataMxd, "Parcels")[0]

point = arcpy.Point(locX, locY)

with arcpy.da.SearchCursor(parcelsLayer, ["OBJECTID","SHAPE@","LEGAL_DESC"] ) as cursor:
    for row in cursor:
        polygonGeom = row[1] # SHAPE@ has index 1 in the fields array above
        if polygonGeom.contains(point):
            print row[0] # Prints OBJECTID value (Field index is 0)
            print row[2] # Prints LEGAL_DESC value (Field index is 2)
            break

[DanPatterson_Retired]

Please mark as answered to facilitate others finding a response that provided useful information.

[IanMurray]

Glad you got it working faster now. 

[TomSellsted]

You might also try this.  If your parcels have a spatial index, this should be pretty fast too.

arcpy.SelectLayerByLocation(in_layer=arcpy.PointGeometry(point), select_features=parcelsLayer)