how can I optimize a script that requires spatial queries, polygon intersects polygon?

LuisEduardoPerezGraterol

Hello

I want to make a script in ArcGis desktop 10.8 for a simple procedure, write for each zone the province that intersects.
The desired result is for the zone layer to have a text field with the names of the province(s) it intersects.

• The zone layer has 1008 records
• the province layer has 51 records.

So far, the options have not been optimal, I detail them below

These are my questions:

1. Is there an optimal way to perform this procedure and to program it?
2. Does the evaluation of the spatial relations for example: polygon.contains(polygon) use the spatial index or should I activate it?
3. How can I perform this process using the spatial index?

For this I am writing two scripts, the first attempt was:
- create a cursor on the zones layer
- select the feature with using its ID with arcpy.SelectLayerByAttribute_management
- then select by location with the provinces
- apply a cursor over the selected to get the list of provinces that will write in the attribute

The problem: this takes too long

Second attempt:
- I create an empty dictionary
- apply a cursor on zones getting the shape field and OBJECTID
- Iterate with a for
- for each case I create a province cursor
- iterate on the provinces
- I evaluate which provinces intersect the zone, in the dictionary I store the id of the zone as a key with a list of the provinces. For this I evaluate if the geometry of zones intersects the province geometry, For example: zona.overlap(province)

The problem: it still takes a long time, 1.8 minutes for 40 zone and they are 1008

JohannesLindner

Cursors do take quite some time, so you should minimize them. You only need to instantiate a cursor 3 times:

1. read the zones
2. read the provinces
3. update the zones

You can do everything else without cursors:

# [ [Zones.OBJECTID, Zones.Shape] ]
zones = [ [oid, shp] for oid, shp in arcpy.da.SearchCursor("Zones", ["OBJECTID", "SHAPE@"])]
# [ [Provinces.Name, Provinces.Shape] ]
provinces = [ [name, shp] for name, shp in arcpy.da.SearchCursor("Provinces", ["Name", "SHAPE@"])]

# {Zones.OBJECTID: Zones.ProvinceNames}
zone_dict = dict()

for z_oid, z_shp in zones:
    p_names = [p_name for p_name, p_shp in provinces if z_shp.overlaps(p_shp)]
    zone_dict[z_oid] = ", ".join(p_names)

with arcpy.da.UpdateCursor("Zones", ["OBJECTID", "ProvinceNames"]) as cursor:
    for oid, names in cursor:
        try:
            names = zone_dict[oid]
            cursor.updateRow([oid, names])
        except KeyError:
            print(f"could not find any province names for zone with OID {oid}")