Arcpy get geometry object at location without selecting or deselecting

843
4
09-20-2021 03:19 PM
GavinPaw
New Contributor

I am writing scripts that take in a layers geometry at a specific location and create a new layer based on it. In order to construct the new geometry I need access to the coordinates of the existing geometry, but cant seem to find an easy way to get it in arcpy. what I have now is some mixture of SelectLayerByLocation to select the objects at a location and then operating on selected objects to get a list of points. Is there a better way to get either a list of points from objects that are at a given location, or geometry objects that represent those objects?

0 Kudos
4 Replies
BlakeTerhune
MVP Regular Contributor

Could you clarify how you determine a feature is at a given location? Check out the Geometry object for a way to get coordinates.

0 Kudos
GavinPaw
New Contributor

I am starting with X and Y coordinates which I use to create a pointGeometry object.

I then use this through arcpy.management.SelectLayerByLocation() to select anything in the targeted layer that exists at the same location as the point.

from there its fairly easy to either get the points of selected objects using arcpy.da.FeatureClassToNumPyArray() with explode_to_points or the geometry objects by getting a list of selected FIDs and using a SearchCursor to get the geometry objects.

Ideally I would be able to get the same without the selecting or deselecting anything and without any strange middle-man steps

0 Kudos
BlakeTerhune
MVP Regular Contributor

A spatial join is the only other thing that comes to mind but I would prefer SelectLayerByLocation(). There's a code sample at the bottom of the Geometry doc page that shows using arcpy.Geometry() as the output for CopyFeatures() and having immediate access to the geometry for everything. You could try that after they're selected instead of the numpy array or search cursor.

0 Kudos
JohannesLindner
MVP Frequent Contributor

If you don't want to deal with selections, you can check if the point intersects the geometry objects of the targeted layer. Depending on the object count, this might be a lot slower.

# construct your query geometry
query_geometry = ...
# optional: buffer
query_geometry = query_geometry.buffer(100)

# read the geometries of the target layer
layer_geometries = [r[0] for r in arcpy.da.SearchCursor(layer, ["SHAPE@"])]

# filter by intersection with query_geometry
intersecting_geometries = [g for g in layer_geometries if not g.disjoint(query_geometry)]

Have a great day!
Johannes