ArcPy Algorithm – Exploding Overlapping to New Non-Overlapping Polygons

algorithmarcpyoverlapping-features

Given multiple polygons that overlap in multiple ways, I would like to export from these features all polygons that don't overlap with others, iteratively.

The product would be a number of features with no overlap that, when summed together, make up the original.

The products could then be used as input to Zonal Statistics, and this would be much faster than iterating Zonal Statistics over each polygon.

I have been trying to code this in ArcPy without success.

Does code to do this already exist?

Best Answer

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.

Map of seven polygons

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

Four-coloring of the polygons

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.

Related Question