Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
    
Do you mean you want to 'flatten' the data into a topologically correct set? –  nagytech Aug 26 '12 at 13:10
    
@Geoist ZonalStats requires polygons that don't overlap. When you have an overlapping collection, the obvious but inefficient solution is to loop over the polys and compute zonal stats one by one. It would be more efficient to select a subset of non-overlapping polys, apply zonalstats to them, and iterate. The question asks how to make such selections efficiently. –  whuber Aug 26 '12 at 21:05
    
whuber - I think that @Geoist is suggesting creating a set of non-overlapping polygons from the intersections of the input polygons. Look at this image - (can't post images in comments?). The input is on the left. The entire region is covered by three polygons, each of which intersects both of the others. The only non-overlapping subsets are the singletons and these do not satisfy Gotanuki's requirement that the union fills the space. I think Geoist is suggesting creating the set of non-intersecting regions on the right which is valid for zonalstats –  Llaves Aug 27 '12 at 1:09
    
I think there is some confusion as to what the final product should be. Could you provide an example? My interpretation is you would like the output to be a selection of polygons that do not overlap--while discarding or dissolving the remaining polygons. Are you working with one or many feature classes? –  Aaron Aug 27 '12 at 1:47
    
@whuber That was what I was suggesting, thanks. –  nagytech Aug 27 '12 at 2:39

4 Answers 4

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.

share|improve this answer
    
I'm not sure if this answers the question, or what the question is, but it is a good answer none the less. –  nagytech Aug 28 '12 at 2:01
    
@Geoist Is there some way I could make the illustration clearer or explain the problem better? –  whuber Aug 29 '12 at 21:18
up vote 4 down vote accepted

The methodology recommended by #whuber inspired me to take a new direction, and here is my arcpy solution, in two functions. The first, called countOverlaps, make two fields, "overlaps" and "ovlpCount" to record for each poly which polys overlapped with it, and how many overlaps occurred. The second function, explodeOverlaps, creates a third field, "expl", which gives a unique integer to each group of non-overlapping polys. The user can then export new fc's based on this field. The process is broken into two functions because I think the countOverlaps tool can prove useful by itself. Please excuse the sloppiness of the code (and the careless naming convention), as it's pretty preliminary, but it works. Also make sure that the "idName" field represents a field with unique IDs (only tested with integer IDs). Thank you whuber for providing me with the framework necessary to approach this problem!

def countOverlaps(fc,idName):
    intersect = arcpy.Intersect_analysis(fc,'intersect')
    findID = arcpy.FindIdentical_management(intersect,"explFindID","Shape")
    arcpy.MakeFeatureLayer_management(intersect,"intlyr")
    arcpy.AddJoin_management("intlyr",arcpy.Describe("intlyr").OIDfieldName,findID,"IN_FID","KEEP_ALL")
    segIDs = {}
    featseqName = "explFindID.FEAT_SEQ"
    idNewName = "intersect."+idName

    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal] = []
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal].append(idVal)

    segIDs2 = {}
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        segIDs2[idVal] = []

    for x,y in segIDs.iteritems():
        for segID in y:
            segIDs2[segID].extend([k for k in y if k != segID])

    for x,y in segIDs2.iteritems():
        segIDs2[x] = list(set(y))

    arcpy.RemoveJoin_management("intlyr",arcpy.Describe(findID).name)

    if 'overlaps' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'overlaps',"TEXT")
    if 'ovlpCount' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'ovlpCount',"SHORT")

    urows = arcpy.UpdateCursor(fc)
    for urow in urows:
        idVal = urow.getValue(idName)
        if segIDs2.get(idVal):
            urow.overlaps = str(segIDs2[idVal]).strip('[]')
            urow.ovlpCount = len(segIDs2[idVal])
        urows.updateRow(urow)

def explodeOverlaps(fc,idName):

    countOverlaps(fc,idName)

    arcpy.AddField_management(fc,'expl',"SHORT")

    urows = arcpy.UpdateCursor(fc,'"overlaps" IS NULL')
    for urow in urows:
        urow.expl = 1
        urows.updateRow(urow)

    i=1
    lyr = arcpy.MakeFeatureLayer_management(fc)
    while int(arcpy.GetCount_management(arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",'"expl" IS NULL')).getOutput(0)) > 0:
        ovList=[]
        urows = arcpy.UpdateCursor(fc,'"expl" IS NULL','','','ovlpCount D')
        for urow in urows:
            ovVal = urow.overlaps
            idVal = urow.getValue(idName)
            intList = ovVal.replace(' ','').split(',')
            for x in intList:
                intList[intList.index(x)] = int(x)
            if idVal not in ovList:
                urow.expl = i
            urows.updateRow(urow)
            ovList.extend(intList)
        i+=1
share|improve this answer
1  
To connect this to my solution: your countOverlaps corresponds to the two lines nbrhoods <- sapply(vertices, neighbors); degrees <- sapply(nbrhoods, length) in my code: degrees is the overlap count. Of course your code is lengthier because it reflects most of the GIS analysis that is taken for granted in my solution: namely, that you first identify which polygons overlap, and that at the end you use the solution to output polygon datasets. It would be a good idea to encapsulate the graph-theoretic calculations, so if you ever find a better coloring algorithm, it would be easy to plug in. –  whuber Aug 28 '12 at 16:03
    
thanks whuber, your solution was a bit more elegant =D –  gotanuki Aug 28 '12 at 17:46

I find it hard to deduce what solution was originally required (before trying to de-overlap intersecting polygons), but had my own bump with zonal stats not allowing overlap. I'll describe my project and solution.

I want to deduce the height of buildings from a very finely meshed elevation grid and building-footprints as polygons.

Since the buildings do not overlap, applying (zonal) stats to deduce elevation-information can be done in a straightforward manner.

However, to determine building-height, I also need the terrain-elevation near the buildings. To extract this information, I buffered the buildings (outside only). These polygons routinely overlap.

Since I really want to apply stats to the gridvalues falling within each polygon, it is not required that the polygons be de-overlapped, it is preferrable that they stay intact (as much as possible), otherwise I would have to reassemble the stats.

Rethinking, I could have written a script selecting the elevations per polygon and processing them to produce stats. I did something similar by writing a C program (standalone executable). It has very good performance, but needs the grid as a .bil and the polygons as shapefile. The stats are similar to zonal stats. Complex and multipart polygons are supported, but no extensive testing has been done sofar.

If anyone thinks this potentially useful, please let me know. Jan

share|improve this answer
1  
Often you can obtain information for terrain near the buildings, without any buffering, by performing a local statistical summary and then summarizing that by zones. –  whuber May 28 '13 at 19:19

In this cases I generally use the following method:

  • Pass the Feature class through a UNION; (It breakes the polygons in all its intersections)
  • Add X, Y and Area fields and calculate them;
  • Dissolve the result by X, Y, Area fields.

I believe the result will be the one you wanted, and you can even count the number of overlaps. Don't know if in terms of performance will it be better for you or not.

share|improve this answer
2  
this method doesn't get you to the desired product, which is a minimal series of selections or unique feature classes of the original that do not overlap. the products will be fed into zonal statistics, and therefore keeping the original geometry of each feature is vital. –  gotanuki Aug 29 '12 at 17:32
    
You are right, sorry. I did not understood the question well. In that case, and depending on the size of the raster, I would normally convert the raster to a temporary point feature class (each cell a point), and perform a spatial join between it and the polygon layer. Maybe its a very simplistic, and performance unfriendly approach but works and the overlapped polygons won't give you any problem. –  Alexandre Neto Aug 29 '12 at 17:58
    
If I understand correctly what you mean by this spatial join, your second solution still won't work, Alexandre, because there is a many-to-many relation between the points and the polygons. Regardless, for any sizable raster this vector-based approach will be extremely inefficient and for large rasters it will be impossible to carry out. –  whuber Aug 29 '12 at 21:16
    
@whuber You are right about being a very slow process (Toke me around half-hour with a 4284 x 3009 raster and 2401 polygons, in a dualcore 2.8Ghz, 3Gb RAM with vista). But It works, as I have tested it already. In the Spatial Join you have to use one to one relation, and aggregate the raster values (as mean, Sum, etc...). The result will be a vector polygon layer similar to the original but with a new column with the aggregated raster values that intersect each polygon. Not being a optimal solution this could be useful for someone with less programming skills (like me :-)). –  Alexandre Neto Aug 30 '12 at 11:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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