After posting my comment yesterday, I realized there was a LOT of room for optimizing the simulations since the polygons and buffer distances do not change during the runs. Since spatial operations are more expensive, often considerably more, than other mathematical operations; identifying all of the spatial interactions first and caching them allows for a spatial analysis to be conducted without using spatial operators. The efficiency of Python sets can be used instead of relying on geoprocessing tools or ArcPy Geometry methods.
I ran the following code against the Redlands, CA zoning data. Instead of picking a region of Redlands, I used the entire city.
# Setup parameters
fc = # path to feature class or shape file
buffer_dist = 500 # buffer distance in units of feature class projection
subset_sql = "CATEGORY = 'Single Family Residential'" # SQL WHERE clause defining subset of units to analyze
num_runs = 100000 # number of scenarios to model
# Buffer and spatial join to identify and cache spatial interactions between
# adjacent/neighboring units
subset = arcpy.MakeFeatureLayer_management(fc, "subset", subset_sql)
subset_buff = arcpy.Buffer_analysis(subset,
"in_memory/buff",
buffer_dist)
subset_rel = arcpy.SpatialJoin_analysis(subset,
subset_buff,
"in_memory/rel",
"JOIN_ONE_TO_MANY",
"KEEP_COMMON")
buff_rel = defaultdict(set)
with arcpy.da.SearchCursor(subset_rel, ["TARGET_FID", "ORIG_FID"]) as cur:
for target, orig in cur:
buff_rel[target].add(orig)
arcpy.Delete_management(subset_buff)
arcpy.Delete_management(subset_rel)
# Retrieve a list of unique/object identifiers
arcpy.SelectLayerByAttribute_management(subset, "SWITCH_SELECTION")
fids = [int(fid)
for fid
in arcpy.Describe(subset).FIDset.split(";")]
arcpy.SelectLayerByAttribute_management(subset, "CLEAR_SELECTION")
# Run simulations/scenarios against the data set
results = set()
for _ in range(num_runs):
shuffle(fids)
sel_set = {fids[0]}
for fid in fids:
if sel_set.isdisjoint(buff_rel[fid]):
sel_set.add(fid)
results.add(frozenset(sel_set))
# Create a size distribution of results
count_hundred = Counter(len(run) for run in results)
Pre-caching the spatial interactions worked better than I expected. I was able to run 100,000 scenarios on the 857 Single Family Residential polygons in around 3 minutes. This approach was several orders of magnitude faster than my original approach.
A breakdown of the results looks like:
Given the number of analysis units, all 100,000 scenarios were unique. Also, there was a single run with 168 polygons identified, the maximum. The maximum is represented by 1 of 100,000 scenarios considered.
