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

LuisEduardoPerezGraterol — Tue, 28 Jun 2022 00:59:04 GMT

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:

Is there an optimal way to perform this procedure and to program it?
Does the evaluation of the spatial relations for example: polygon.contains(polygon) use the spatial index or should I activate it?
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

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

JohannesLindner — Tue, 28 Jun 2022 10:23:17 GMT

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

read the zones
read the provinces
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}")

topic Re: how can I optimize a script that requires spatial queries, polygon intersects polygon? in ArcGIS Pro Questions

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

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