Using ArcGIS Desktop 10.6, I have the following workflow in a python script to try and calculate coverage area affected on a wireless network if certain towers went down:

  1. User enters a list of polygon IDs (these are the tower outages)
  2. Using an in-memory layer, script selects all polygons from a layer in a file geodatabase, except the ones listed by the user

Important note: Polygon layer is 500GB in size and cannot be generalized.
Contains complex geometries and often intentional overlap between polygons

  1. Polygons from #2 are merged (arcpy.Merge_management)

  2. Loop through all the polygons that were removed and clip any overlapping areas from merged result in #3

  3. Merge output from #4 and calculate its area.

I am fairly confident that I have optimized the script as much as possible, and due to the size of the data, it is taking a significantly large time to process a result.

The question I have is what is a good alternative approach that will significantly reduce the processing time for the user but still accomplish the same workflow?

Things I am currently contemplating:

  • Pre-cooking every combination of result using a union on the original layer. This may be a beast to process, but would only need to happen every ~6mths.
  • Storing in SQL Geometry and using native SQL spatial functions to return results. Not tested, but I assume some of the more complex geoms might cause some issues.
  • Raster based approach? Not put a lot of thought into this one.
  • When you loop through all the removed polygons are you looping through all or selecting by location and looping through the selection... no point clipping something that doesn't overlap. If you want performance you might need to code this in C# or VB.net for ArcObjects.. it's a simple workflow but would be a lot faster in ArcObjects. – Michael Stimson May 30 '18 at 2:05
  • I'm thinking that a raster approach would be worth looking at. Replace 3. with a PolygonToRaster_conversion. Do the same with the user list of polygon IDs, then subtract one from the other and convert the result back into polygon (if required). You'll lose some accuracy in the area calculation (depending on cellsize) but that may be acceptable if the processing time is significantly reduced. – Dan May 30 '18 at 2:14
  • How are you performing the select operation? – Aaron May 30 '18 at 3:24
  • I will try and cut this down a bit later, but here is the current script in its raw form pastebin.com/9Q4Cuc4H l.283 is the function in question – Simon May 30 '18 at 3:38
  • 1
    @BERA that is an approach I think I will test. Never tried a union on a 500GB layer, as long as it works, the disk space should not be an issue. I also like the raster approach, and further upstream the data actually originates in raster, so will also be investigating that. I think I have oversimplified the workflow a bit. The linked script has extra info for dealing with a variety of layers (3G/4G/LTE) and a whole other scenario. Will attempt to simplify this down. – Simon May 30 '18 at 8:37

Not sure I understand you completely, you want to find offline areas, not covered by any of the online antennas? Try the Union alternative you are thinking about and then Select By Attributes. The Union will be pretty slow but not as slow as merge-clip-merge every time you want to run the analysis and you only need to Union every six months. And SQL queries are fast.


Split by Attributes (Object ID or some other ID) and then Union all outputs: enter image description here Construct a Query using Python, select and sum:

import arcpy

fc = r'C:\Arenden\Default.gdb\Network_areas'
splitfield = 'SplitID'

offline_towers = [10,11,12]

offline_fields = ['FID_T'+str(f) for f in offline_towers]
online_fields = ['FID_T'+str(f[0]) for f in arcpy.da.SearchCursor(fc,splitfield) if f[0] not in offline_towers]

sql_list = []
for field in offline_fields:
    sql_list.append("{0}>-1 AND {1}".format(field,
                                ' AND '.join(['{0}<0'.format(f) for f in online_fields])))

sql = ' OR '.join(sql_list)

print sum([i[0] for i in arcpy.da.SearchCursor("Network_areas_Union",'SHAPE@AREA')])

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.