I am facing a similar problem to Using Zonal Statistics As Table for overlapping polygons in ArcPy? which links to Calculating zonal statistics of raster data in multiple overlapping zones and combining them into one table and an Esri supplementary toolbox.
I have a number of building polygons which need to have the statistical information from a raster attributed for a buffer zone and these buffers frequently overlap. All these previous questions' solution include iterating the geometries and generating a zonal statistic for each individual input polygon; I did have some hope that the source code for the ZonalStatisticsAsTable2 from Esri would have a different solution as it appeared to planarize the polygons but then went on to perform a zonal statistic for each fragment, lines 289-309:
# Perform zonal statistics for each class
temp_lyr = "temp_layer"
cl_separator = ' OR \"%s\" = ' % oid_field
for index, cl in enumerate(classes):
arcpy.SetProgressorLabel(
"Processing layer %d of %d..." % (index+1, num_classes))
where_clause = '\"%s\" = %s' % (oid_field, \
cl_separator.join(map(str, classes[cl])))
temp_table = os.path.join(temp_dir, "zone_%d.dbf" % index)
arcpy.MakeFeatureLayer_management(temp_features, temp_lyr, \
where_clause)
try:
arcpy.sa.ZonalStatisticsAsTable(temp_lyr, feature_field, \
value_file, temp_table, ignore_value, statistic)
except:
arcpy.GetMessages(0)
# Merge tables
arcpy.env.workspace = temp_dir
table_list = arcpy.ListTables("zone*")
arcpy.Merge_management(table_list, output_table)
del table_list
The root cause of the problem comes from the method used by Zonal Statistics as Table
If the zone input is a feature dataset, a vector-to-raster conversion will be internally applied to it. To ensure that the results of the conversion will align properly with the value raster, it is recommended that you check that the extent and snap raster are set appropriately in the environment settings and the raster settings.
Overlapping polygons overwrite with the zone id of the latter polygon, in some cases obliterating completely, in others the buffers are truncated and no longer accurate. Assuming that iterating each input is at least a workable solution I have tested on a small dataset of about 2k buildings and find it forbiddingly slow, far to slow to roll out:
ZonalDict={}
with arcpy.da.SearchCursor(BuffA_FC,'SourceID') as bCur:
for bRow in bCur:
LayerID = 'Feat_{}'.format(bRow[0])
ZS_Name = 'IN_MEMORY\\Tab_{}'.format(bRow[0])
LayerDQ = 'SourceID = {}'.format(bRow[0])
BuffALyr = arcpy.MakeFeatureLayer_management(BuffA_FC,LayerID,LayerDQ) # make feature layer should be quicker than exporting features with Select
arcpy.sa.ZonalStatisticsAsTable (BuffALyr,'SourceID',ShrubRasterPcnt,ZS_Name,statistics_type='MEAN')
with arcpy.da.SearchCursor(ZS_Name,'MEAN') as zsCur:
for zsRow in zsCur:
ZonalDict[bRow[0]]=zsRow[0] # embed the mean for this feature in the dict with key of source identifier
# cleanup: important if arcpy.env.overwriteOutput is not set to True
arcpy.Delete_management(ZS_Name)
arcpy.Delete_management(LayerID)
with arcpy.da.UpdateCursor(BuildingOutlines,['OID@',OutputFields[Shrub_A_Field]]) as UCur:
for uRow in UCur:
uRow[1] = ZonalDict[uRow[0]] # I should not need to implement if uRow[0] in ZonalDict, it should be guaranteed
UCur.updateRow(uRow)
I have an idea that a faster option should exist, perhaps planarizing the polygons and creating a lookup from the new IDs of the planar polygons to the source overlapping polygons then a simple Zonal Statistics as Table should be sufficient with a double join with a statistic of the fragment statistics but I'm having difficulty generating substance from this nebulous idea which may mean it's unworkable. Or perhaps identifying the overlapping buffers, segregating and iterating the overlapping buffers then appending to the statistical table generated directly from the disparate buffer data.. gut feel is that identifying the overlapping polygons by iteration could be just as slow.
Does anyone have an idea that overcomes the overlapping zone problem that isn't so slow as to be unworkable?