Iterate selection and calculate field not working

[ThomasColson]

With the following, I am trying to select each polygon by its name, select points from another feature class that are within the selected polygon, then calculate a field using the same name that was used to select the polygon (by attributes). A spatial intersect, or any other method that results in creating a new FC as output and having to do joins won't work here, as this is (hopefully) to be an unattended nightly update script. It iterates through the polygons fine, but the result is that ALL points get the name of the last polygon. 

import os
import errno
import arcpy
import sys  
import traceback  
from arcpy import env
update_feature_class = r'C:\Users\tcolson\Documents\ArcGIS\Projects\MyProject\MyProject.gdb\GRSM_RESEARCH_LOCATIONS'
county_feature_class = r'C:\Users\tcolson\Documents\ArcGIS\Projects\MyProject\MyProject.gdb\CountyBoundaries'

with arcpy.da.SearchCursor(county_feature_class,['SHAPE@','NAME']) as cursor:
    for row in cursor:
        print('Selecting '+row[1])
        expression = "NAME = '"+row[1]+"'"
        arcpy.SelectLayerByAttribute_management (county_feature_class, "NEW_SELECTION", expression)
        arcpy.SelectLayerByLocation_management(update_feature_class, "INTERSECT", county_feature_class)
        arcpy.CalculateField_management(update_feature_class, "COUNTY","'"+row[1]+"'", "PYTHON3")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[ThomasColson]

Here is something interesting: when I try to extend the point-in-polygon script to polygon-centroid-in-polygon (what administrative boundary is this building located in?), if the input polygon intersects the within polygon (building straddles a county boundary), the result is null: 

    with arcpy.da.UpdateCursor(update_feature_class, ['SHAPE@XY', update_name]) as update_cursor:

        # Iterate through points
        for point in update_cursor:
            point_count += 1
            name_list = []
            for name in poly_names:
                comparison_count += 1
                poly = relateDict[name]
                # Check if the current point is within the current polygon
                if point[0].within(poly):  # Returns True or False
                    # Add the polygon name to the list
                    name_list.append(name)
            point[1] = ",".join(name_list)
            # Update the cursor with the new value
            update_cursor.updateRow(point)
    arcpy.AddMessage("Processed {0} Polygon/Polygon or Line comparisons for {1} Polygons and {2} Polygons or Lines".format(comparison_count, poly_count, point_count))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[JoshuaBixby]

SHAPE@XY doesn't return a point geometry, or even a point, it returns a tuple containing coordinates.  In your code snippet as posted, the point[0].within(poly) line should generate an AttributeError since tuples don't contain a within method.  I am guessing the code snippet is abridged and something got left out when you pasted it; otherwise, I don't understand how the code is running at all.

[ThomasColson]

Hmmmm   How To: Calculate feature centroids  suggests:

"In the following sample, the Data Access Module is used to retrieve the centroid coordinates using the SHAPE@XY token." which I interpret as the polygon centroid. If that is indeed a tuple of many coordinates, how does one get the centroid using the da method?

[DanPatterson_Retired]

Polygon—ArcPy classes | ArcGIS Desktop 

using SHAPE@ to get the polygon

then

centroid(Read Only)The true centroid if it is within or on the feature;
 otherwise, the label point is returned. Returns a point object. ‍‍‍‍‍

so

cent = row[whatever].centroid # which returns a point object

cent.X and cent.Y if you need the coordinates of the centroid and not the point object itself

[JoshuaBixby]

I suggest using Dan's suggestion and returning the geometry itself and then calling the centroid method.

Regarding documentation, the next sentence after the one you quoted is, "These coordinates are then used to build point geometry,...."  As helpful as Esri Support can be, in the hierarchy of documentation they are a step or two below the documentation from Esri Development.  Looking at UpdateCursor—Help | ArcGIS Desktop , the object type is clearly documented:

• SHAPE@XY —A tuple of the feature's centroid x,y coordinates.
• SHAPE@TRUECENTROID —A tuple of the feature's true centroid x,y coordinates.

[ThomasColson]

Here's the full snippet (which I've attached to a toolbox). Even using SHAPE@ I still get the nulls for polygons that straddle a polygon boundary. Have also tried  SHAPE@ and SHAPE@TRUECENTROID

import arcpy
update_feature_class = arcpy.GetParameterAsText(0)
poly_feature_class = arcpy.GetParameterAsText(1)
update_name = arcpy.GetParameterAsText(2)
poly_name = arcpy.GetParameterAsText(3)

poly_count = 0
relateDict = {}
# Set up search cursor for poly fc
with arcpy.da.SearchCursor(poly_feature_class,[poly_name,'SHAPE@']) as search_cursor:

    # Iterate through polygons
    for polygon in search_cursor:
        if not polygon[0] in relateDict:
            relateDict[polygon[0]] = polygon[1]
        else:
            relateDict[polygon[0]] = relateDict[polygon[0]].union(polygon[1])

poly_names = []
for name in sorted(relateDict):
    poly_names.append(name)
    poly_count += 1

point_count = 0
comparison_count = 0
# Set up update cursor for point fc
with arcpy.da.UpdateCursor(update_feature_class, ['SHAPE@', update_name]) as update_cursor:

    # Iterate through points
    for point in update_cursor:
        point_count += 1
        name_list = []
        for name in poly_names:
            comparison_count += 1
            poly = relateDict[name]
            # Check if the current point is within the current polygon
            if point[0].within(poly):  # Returns True or False
                # Add the polygon name to the list
                name_list.append(name)
        point[1] = ",".join(name_list)
        # Update the cursor with the new value
        update_cursor.updateRow(point)
arcpy.AddMessage("Processed {0} Polygon/Point comparisons for {1} Polygons and {2} Points".format(comparison_count, poly_count, point_count))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

How do they straddle a boundary? is that determined via polygon to polygon intersection? (ie polygon intersects polygon …. not polygon within polygon)

If you are comparing the centroid of a straddling polygon, the centroid may not be within the bounds of the intersecting polygon

[ThomasColson]

[DanPatterson_Retired]

centroid within a distance of.... if you are doing some kind of select by attributes (thread is too long and convoluted now)

[JoshuaBixby]

I agree with Dan that either starting a new question or branching this question to another thread would be good given the length and slight topic change.

What I will say here, though, is that reading through the script for me is made more difficult by naming rows from the cursors with names like "polygon" and "point."  Typically when something is named "polygon," it is a polygon object and not a row containing a polygon object and other attributes.  Maybe it is just me, but I catch myself having to go back and remind myself what exactly "polygon" and "point" are when I see polygon[0] because ArcPy polygon objects do support indexing and slicing.

Beyond variable naming, I am trying to understand why you are creating a poly_count counter when you can just check the length of the poly_names list to get the count.  And I don't see poly_count actually being used anywhere.

Since Python dictionaries are not inherently spatial, looping over them and doing spatial operations on each geometry stored within them can be slow depending on how many polygons are being iterated over, but I haven't followed this thread closely until now, so maybe that discussion already happened.  If so, disregard this last comment.