This is a graph coloring problem.
Recall that a graph coloring is an assignment of a color to the vertices of a graph in such a way that no two vertices which share an edge will also have the same color. Specifically, the (abstract) vertices of the graph are the polygons. Two vertices are connected with an (undirected) edge whenever they intersect (as polygons). If we take any solution to the problem--which is a sequence of (say k) disjoint collections of the polygons--and assign a unique color to each collection in the sequence, then we will have obtained a k-coloring of the graph. It is desirable to find a small k.
This problem is pretty hard and remains unsolved for arbitrary graphs. Consider an approximate solution that's simple to code. A sequential algorithm ought to do. The Welsh-Powell algorithm is a greedy solution based on a descending ordering of the vertices by degree. Translated to the language of the original polygons, first sort the polygons in descending order of the number of other polygons they overlap. Working in order, give the first polygon an initial color. In each successive step, try to color the next polygon with an existing color: that is, choose a color that is not already used by any of that polygon's neighbors. (There are many ways to choose among the available colors; try either the one that has been least used or else choose one randomly.) If the next polygon cannot be colored with an existing color, create a new color and color it with that.
Once you have achieved a coloring with a small number of colors, perform zonalstats color by color: by construction, you're guaranteed that no two polygons of a given color overlap.
Here's sample code in R. (Python code wouldn't be much different.) First, we describe overlaps among the seven polygons shown.
edges <- matrix(c(1,2, 2,3, 3,4, 4,5, 5,1, 2,6, 4,6, 4,7, 5,7, 1,7), ncol=2, byrow=TRUE)
That is, polygons 1 and 2 overlap, and so do polygons 2 and 3, 3 and 4, ..., 1 and 7.
Sort the vertices by descending degree:
vertices <- unique(as.vector(edges))
neighbors <- function(i) union(edges[edges[, 1]==i,2], edges[edges[, 2]==i,1])
nbrhoods <- sapply(vertices, neighbors)
degrees <- sapply(nbrhoods, length)
v <- vertices[rev(order(degrees))]
A (crude) sequential coloring algorithm uses the earliest available color not already used by any overlapping polygon:
color <- function(i) {
n <- neighbors(i)
candidate <- min(setdiff(1:color.next, colors[n]))
if (candidate==color.next) color.next <<- color.next+1
colors[i] <<- candidate
}
Initialize the data structures (colors and color.next) and apply the algorithm:
colors <- rep(0, length(vertices))
color.next <- 1
temp <- sapply(v, color)
Split the polygons into groups according to color:
split(vertices, colors)
The output in this example uses four colors:
$`1`
[1] 2 4
$`2`
[1] 3 6 7
$`3`
[1] 5
$`4`
[1] 1
It has partitioned the polygons into four non-overlapping groups. In this case the solution is not optimal ({{3,6,5}, {2,4}, {1,7}} is a three-coloring for this graph). In general the solution it gets shouldn't be too bad, though.
