How to obtain the value of a feature programatically, based on a coordinate

Question

I have a point layer, and a polygon layer, and I wish to calculate the value of an attribute from the polygon layer at each point's location. Eg, find the area of the polygon which is found at each point, and write it to the point's attribute table:

One simple approach is to perform a spatial join - but this creates a new output point layer. I would prefer to update the existing point layer.

I could iterate through the points layer and:

1. obtain each point's OBJECTID
2. create a layer from the selected point, using the OID as the Where clause
3. perform a Select By Location, using the above layer as the select_features
4. iterate through the selected polygon(s) and write the Area value to the point

Steps 2 and 3 seem a bit extraneous - is there a programmatic way to obtain the value of a layer based solely on an XY coordinate?

This would be akin to using the Identify tool in ArcMap, whereby clicking on a location reveals the attributes of the layers at that location.

Answer 1

You can just iterate over the points with an UpdateCursor, make a geometry object out of the SHAPE@ and then perform steps 3 and 4:

  with arcpy.da.UpdateCursor(points, ["SHAPE@", "AREA_FROM_POLYGON"]) as cur:
        for row in cur:
            xy = row[0]
            arcpy.SelectLayerByLocation_management(polygons_lyr, "CONTAINS", xy)
            if arcpy.da.SearchCursor(polygons_lyr, ["AREA"]):
                with arcpy.da.SearchCursor(polygons_lyr, "[AREA"]) as cur_polyg:
                    for row_polyg in cur_polyg:   ## this assumes there's only 1 polygon intersecting a point
                        area = row_polyg[0]  
                row[1] = area
                cur.updateRow(row)
            else:
                arcpy.AddMessage("Point doesn't intersect any polygon")

I don't think there's a more straightforward way if you don't want to copy the points layer. If your polygon layer was a raster, you could use Get Cell Value with the SHAPE@XY, but this is maybe out of scope for you.

Answer 2

One possible workaround is to convert your polygon layer to raster based on the area field, and use "extract multiple value to point". This will create an intermediate file(the raster) and would be an issue for a large number of fields, but you don't create a new point feature.

Another workaround is to join the points to the polygons, then join the attribute table of the resulting layer using the OID of the points. Again, this create an intermediate file but you bon't duplicate your point feature.

I don't see other methods without looping on each point.