Writing polygon values to point feature class based on an intersection of geometry.

Question

Hello, I want to create an automated scheduled task with Notebooks/Python working with hosted feature layers, where polygon values are written to point features based on an intersection of geometry. I have 2 point feature classes, let's call them: PointFC1 PointFC2 and 2 polygon feature classes, let's call them: PolyFC1 PolyFC2 I would like to write PolyFC1 'Location' field to PointFC1 and PointFC2 'Location' field, and write PolyFC2 'Zone' to PointFC1 and PointFC2 'Zone' fields based on whether those point features fall within the boundary of said polygons. I would only like to run this calculation on points where 'Location' or 'Zone' is '', so redundant processing is removed. I have done this in ArcGIS Pro with field calculations via Python on the individual layers with: getZone(!Shape!) ------------------------------------------------------ import arcpy def getZone(point_geom): polygon_layer = "polygon_layer" # just the layer name as it appears in ArcGIS Pro with arcpy.da.SearchCursor(polygon_layer, ["SHAPE@", "field_name"]) as cursor: for row in cursor: if row[0].contains(point_geom): return row[1] return None I now want to group all the field calculations and run an automated task on them fortnightly, while only doing the points that are in those fields, as both point feature classes link to a Survey123 form where data is collected daily. That way, only features that need to be processed are processed. Any help on this would be much appreciated. Cheers.

Hayden Welch · Answer

I prefer to use Cursors for all my joins and relates since in my experience they're much faster than any other option. SearchCursors are cheap and you can use them to extract and manipulate data all you want, then structure that and apply it in one UpdateCursor operation:

from collections.abc import Iterator, Sequence
from typing import Any

from arcpy import Polygon
from arcpy.da import SearchCursor, UpdateCursor

type UpdateMap = dict[int, dict[str, Any]]
"""Mapping of OBJECTID to Field Updates {FieldName: FieldValue, ...}"""

type FeatureMap = dict[str, str]
"""Mapping of FeatureClass to extracted field"""

def get_shapes(fc: str, field: str) -> Iterator[tuple[Polygon, str]]:
    yield from SearchCursor(fc, ['SHAPE@', field])

def update_feats(fc: str, updates: UpdateMap, fields: Sequence[str]) -> None:
    """Applies an UpdateMap to the target FeatureClass (Only supplied fields are updated)"""
    with UpdateCursor(fc, ['OID@'] + list(fields)) as cur:
        for oid, _ in filter(lambda r: r[0] in updates, cur):
            cur.updateRow([updates[oid][f] for f in cur.fields[1:]])

def get_intersecting(fc: str, field: str, shape: Polygon, value: Any) -> UpdateMap:
    """Gets all features in the input FeatureClass and matches them with 
    the target value based on spatial relation
    """
    return {
        oid: {field: value}
        for oid, 
        in SearchCursor(
            fc, ['OID@'], 
            spatial_filter=shape, 
            where_clause=f'{field} != {value} OR {field} IS NULL',
            search_order='SPATIALFIRST',
        )
    }

def merge_updates(*updates: UpdateMap) -> UpdateMap:
    """Merges multiple UpdateMaps into one"""
    merged: UpdateMap = {}
    for update in updates:
        for oid, vals in update.items():
            if oid in merged:
                merged[oid].update(vals)
            else:
                merged[oid] = vals
    return merged

def apply_join(fc: str, feat_map: FeatureMap):
    """Applies Updates to the input features based on the mapping supplied"""
    updates: list[UpdateMap] = []
    for mapped_feat, field in feat_map.items():
        for shape, val in get_shapes(mapped_feat, field):
            updates.append(get_intersecting(fc, field, shape, val))
    update_feats(fc, merge_updates(*updates), list(feat_map.values()))

def main():
    poly_features: FeatureMap = {
        'PolyFC1': 'Zone',
        'PolyFC2': 'Location',
    }
    apply_join('PointFC1', poly_features)
    apply_join('PointFC2', poly_features)

That apply_join function is the entire business logic, everything else is just to abstract the cursor operations away.

1) Iterate the PolygonClass/Value mapping (PolyFC1 has Zone, etc.)

2) Iterate all Shape/Value combos in that featureclass

3) Get all intersecting points for each shape and add them to an update mapping

4) Merge all updates into one Update map for the Target Point features

5) Apply the updates

Additionally, if you aren't sure that your target features have a matching field name, you can ensure it with a guard clause in the update_feats function:

def update_feats(fc: str, updates: UpdateMap, fields: Sequence[str]) -> None:
    """Applies an UpdateMap to the target FeatureClass (Only supplied fields are updated)"""

for field in fields:
        ensure_field(fc, field)
    get_fields.cache_clear()

with UpdateCursor(fc, ['OID@'] + list(fields)) as cur:
        for oid, _ in filter(lambda r: r[0] in updates, cur):
            cur.updateRow([updates[oid][f] for f in cur.fields[1:]])

@cache
def get_fields(fc: str) -> list[str]:
    return [f.name for f in ListFields(fc)]

def ensure_field(fc: str, field: str, field_info: dict[str, Any] | None = None, default_type: str = 'TEXT'):    
    if not field_info:
        field_info = {}
    if field not in get_fields(fc):
        AddField(fc, field, field_info.get(field, default_type))

the @cache decorator is from functools and will prevent the ListFields function from running for every input field. We need to clear the cache after ensuring fields since we've now update the field list. If you want, you can use a global map for the field_info parameter, or you can expose that in update_feats and pass it along.

This is definitely a overengineered solution to this, but I always prefer to be able to adjust the main loop when things change without manually adding tons more code blocks.

Hayden Welch · Answer

It's kinda dangerous to nest cursors like that. If something goes wrong, you can end up with borked data.

DavidSolari · Answer

You can do all of this in Python, but working in the Spatial Join tool should be more performant as your datasets grow. For each point FC you can join in both polygon FCs, then loop over the results table in a Search Cursor to build up a dictionary of TARGET_FID to Location & Zone, then plug the Object ID of each point into that dictionary during the Update Cursor. This is what my Spatial Join parameters roughly looked like:

arcpy.analysis.SpatialJoin(
    target_features="PT1",
    join_features="PG1 JC1;PG2 JC2",
    out_feature_class=r"memory\PT1_SpatialJoin",
    join_operation="JOIN_ONE_TO_ONE",
    join_type="KEEP_ALL",
    field_mapping='location "location" true true false 255 Text 0 0,First,#,PG1,location,0,254;zone "zone" true true false 255 Text 0 0,First,#,PG2,zone,0,254',
    match_option="INTERSECT"
)