A quicker way to perform zonal statistics for overlapping zone data

Question

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?

Answer 1

The answer came to me overnight, I hadn't hit the 'go' button and seeing as I'd already written the question I might as well finish it off. The process goes like this:

• Planarize the polygons, this can be done with the Esri tool Feature to Polygon, make sure a unique field exists in the planar polygons that is integer and populated with a unique value - Don't rely on FID values.
• Generate centroids for your planar polygons, these will be used later.
• Set your environments of cell size and snap raster to match your data source raster, this speeds up the raster processing as the cells will align which avoids resampling.
• Generate a zone raster from the planar polygons, for very large datasets with more than 2,147,483,647 polygons after planarizing some tiling or batching might be needed.
• Generate a zonal statistics as table using the zone raster and the data source raster.
• Join the statistics to the centroids with Join Field which is permanent.
• Overlay your original buffers with a spatial join or intersect, this will duplicate any centroids that appear in multiple buffers.
• The tricky bit is to generate the required statistic from multiple rows, I used Summary Statistics on the joined centroids to get the sum of sum (from zonal statistics) and sum of count (from zonal statistics) with a case field of the planar polygon unique identifier. The statistic I was after is mean which is sum of sum divided by sum of count, minimum and maximum should also be able to be calculated this way but majority and standard deviation can't.
• Join the summary table to the original buffers... that's it, all done!

The concept of this method could be utilized by any software package if the required tools can be identified; I am using ArcGIS Desktop and so have used arcpy, for the benefit of arcpy users, who have both an advanced license seat and Spatial Analyst extension license the code is:

arcpy.env.snapRaster = DataRaster
arcpy.env.extent     = DataRaster
arcpy.env.cellSize   = DataRaster

WorkDB    = 'IN_MEMORY' # Can be changed to a GDB if persisting for debugging is required
BuffA_FC  = os.path.join(WorkDB,'BldBuf_A') 
A_Pol     = os.path.join(WorkDB,'A_Pol')
A_Pol_Ras = os.path.join(WorkDB,'A_Pol_Ras')
A_Point   = os.path.join(WorkDB,'A_Point')
A_Over    = os.path.join(WorkDB,'A_Point_Poly_Overlay')
A_OverTab = os.path.join(WorkDB,'A_OverlayStatistics')
A_ZSTab   = os.path.join(WorkDB,'ZS_BuffA')

arcpy.AddMessage('Part A {} metres'.format(BuffADist))
# Prior to this the building features have been attributed with a field 'SourceID' which is
# a copy of the arcpy.Describe(BuildingFC).OIDFieldName field to tie all the layers together
arcpy.FeatureToPolygon_management(BuffA_FC,A_Pol,attributes='NO_ATTRIBUTES') # Planarize polygons

# Planarize the polygons to resolve overlapping areas as their own polygon. The variable A_Pol is
# the buffer polygon feature class generated from BuildingFC feature class
arcpy.PolygonToRaster_conversion(A_Pol,arcpy.Describe(A_Pol).OIDFieldName,A_Pol_Ras,'MAXIMUM_AREA')

# Two step process for the tighter control of the joined fields, too many fields makes debugging a pain
# create a value raster using the FID of the planar polygons, the values don't matter at this point
# but must be integer and unique which fits FID, OID or OBJECTID fields
arcpy.sa.ZonalStatisticsAsTable (A_Pol_Ras,'VALUE',DataRaster,A_ZSTab)
arcpy.JoinField_management(A_Pol,arcpy.Describe(A_Pol).OIDFieldName,A_ZSTab,'VALUE',['SUM','COUNT'])

# Generate points and spatial join to the buffer polygons. Point centroids are use to avoid sliver areas 
# and false joins then use summary statistics to compile the sum and count of cells with a case unique
# field from the overlay to associate the summary table to the buffers, if a point falls within two
# buffers it will be duplicated in the output points allowing the associated area to be attributed
# to all overlapping polygon buffers.
arcpy.FeatureToPoint_management(A_Pol,A_Point,'INSIDE')
arcpy.SpatialJoin_analysis(A_Point,BuffA_FC,A_Over,'JOIN_ONE_TO_MANY','KEEP_COMMON')
arcpy.Statistics_analysis(A_Over,A_OverTab,[['COUNT','SUM'],['SUM','SUM']],'SourceID' )
arcpy.AddField_management(A_OverTab,'Mean_Cover','DOUBLE')

# calculate the statistic for each unique polygon now that the total sum and count for
# each input buffer polygon has been ascertained
arcpy.CalculateField_management(A_OverTab,'Mean_Cover','!SUM_SUM! / !SUM_COUNT!',"PYTHON")

# If all you are interested in is getting the value for the buffer stop here
arcpy.JoinField_management(BuffA_FC,'SourceID',A_OverTab,'SourceID',"Mean_Cover")

# My goal is to attribute the source building with the value of the buffer so extra steps are needed
# Populate a dictionary with the source ids and mean cover to identify with source building polygons
ZonalDict={}
with arcpy.da.SearchCursor(A_OverTab,['SourceID','Mean_Cover']) as bCur:
    for bRow in bCur:
        ZonalDict[bRow[0]]=bRow[1]

# Update the mean_cover field with the values in the dictionary being sure to confirm
# that the source ID exists. If there is no matching ID in the dict either the tool
# is being run with a subset of source feature or something is horribly wrong
with arcpy.da.UpdateCursor(BuildingFC,['OID@','Buff_Mean') as UCur:
    for uRow in UCur:
        if uRow[0] in ZonalDict:
            uRow[1] = ZonalDict[uRow[0]]
        else:
            uRow[1] = -1 # Error flag value, indicating no matching value in the dict for this key
        UCur.updateRow(uRow) # don't forget to store the row.

Answer 2

This question is repeated so many times in the past and I am one of those who asked this question. After looking for a solution for this problem, my definite answer is as follows:

1. QGIS can calculate the zonal statistics with overlapping polygon.

2. R can calculate the zonal statistics. For example. exact_extract(raster_object,polygon_object c('mean'))) can do what you want.

I hope that those tips are helpful.

Answer 3

Apparently now the same tool of the ArcGIS Pro version Zonal Statistics As Table does iterate each part when there are overlaps.

Notice the differences in the tool description for ArcMap (v10.8):

If the Input raster or feature zone data (in_zone_data in Python) has overlapping polygons, the zonal analysis will not be performed for each individual polygon. Since the feature input is converted to a raster, each location can only have one value.

An alternative method is to process the zonal operation iteratively for each of the polygon zones and collate the results.

And the tool description of the ArcGIS Pro version (v3.0):

If the Input Raster or Feature Zone Data has overlapping features, the zonal analysis will be performed for each individual feature.