R – How to Dissolve Only Overlapping Polygons Using R

dissolvegeoprocessingmaptoolsrrgeos

I have a spatialpolygons-dataframe in R that contains overlapping polygon features. I would like to dissolve only the overlapping features into separate polygon features (one feature per par of overlapping polygons) and preserve the other features as they are. The length of my dataset should therefore reduce e.g. from 4000 to 3900 if say 200 polygons overlap pairwise.

A simplified visualization:
enter image description here

I tried the unionSpatialPolygons() from the maptools package as indicated here. However this will dissolve all polygons into a single polygon.

A workaround would be to use an id-vector that defines which polygons should be dissolved as indicated by the help manual.

if the id argument is used, it should be a character vector defining the memberships of the output Polygons objects, equal in length to the length of the polygons slot of spgeom

In my case it would be those polygons that overlap together (can be two or more).

Edit

A possible but not efficient solution is given in the code below, assuming you already have an intersection matrix, that you can create e.g. with gIntersects setting byid=TRUE. IF anyone knows an easier and more efficient solution for this, please contribute. 

### example data (this should later be your intersection matrix)

mx<-matrix(c(TRUE,TRUE,FALSE,FALSE,FALSE,FALSE,TRUE,
             TRUE,TRUE,TRUE,FALSE,FALSE,FALSE,FALSE,
             FALSE,TRUE,TRUE,FALSE,FALSE,FALSE,FALSE,
             FALSE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE,
             FALSE,FALSE,FALSE,TRUE,TRUE,FALSE,FALSE,
             FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,FALSE,
             TRUE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE),7)

### groupings
# create a list for the results
results.list<-as.list(vector(length=ncol(mx)))

# group
for(i in 1:ncol(mx)) {
  tmp <- which(mx[,i]) # get TRUE FALSE values for the i-th column
  ifelse(length(tmp)>1, # if there is more than 1 TRUE value,
         tmp.expand<-which(apply(mx[,tmp],1,any)), # get the row-number of the other TRUE Values and create a vector called expand
         tmp.expand<-tmp) # otherwise define tmp as expand
  while(length(tmp.expand)>length(tmp)){ # while tmp.expand has more items than tmp
    tmp<-tmp.expand # reset tmp to tmp.expand
    tmp.expand<-as.vector(which(apply(mx[,tmp],1,any))) # get all new row-numbers of new TRUE values
  }

  results.list[[i]]<-tmp.expand # store results in the list
  print(paste("nr", i, "out of", ncol(mx),"done", sep=" "))
}

# create unique ids from the results
results.list<-
llply(.data = results.list,
      .fun = function(x)paste(x,collapse=""))

Now you can use this list as an ID vector in the unionSpatialPolygons function. This will create new polygons from the overlapping ones and leave the others as they are.

Note
The amount of computational power required by this approach increases exponentially with the size of your matrix/nr of polygons and the nr. of overlappings. If you have a very big data-set you might rather subset it first and than process the subsets separately. After you can join them again and apply the function again to get the same result.
I also tried the code with the lapply function instead of a for loop but it is not really faster, at least if applied on a single core.

The code was partly developed with help of this question on SO.

Best Answer

Here is a short example. I assume that by overlay you are looking for intersects; namely to dissolve all polygons that intersect from both layers.

library(sp)
library(rgeos)
# Create a dataset
poly <- SpatialPolygons(list(
  Polygons(list(Polygon(coords = matrix(c(1, 1, 4, 3, 4, 2, 1, 1), ncol =  2, byrow = TRUE))), ID = "1"),
  Polygons(list(Polygon(coords = matrix(c(3, 1.5, 3.5, 2, 4, 2, 3, 1.5), ncol = 2, byrow = TRUE))), ID = "2"),
  Polygons(list(Polygon(coords = matrix(c(4, 1, 5, 2, 5, 1, 4, 1), ncol = 2, byrow = TRUE))), ID = "3")
))

# Split dataset to two SpatialPolygon objects
polyA <- poly[1, ]
polyB <- poly[2:3, ]

# Show dataset
plot(poly, axes = TRUE)
plot(polyA, axes = TRUE, add = TRUE)
plot(polyB, axes = TRUE, add =TRUE, col = "red")

You can see the sample data set below. PolyB has two polygons from which one intersects the polygon in PolyA.

Dataset basic plot

Using gUnion will result in one polygon of all polygons in both objects, as you have suggested: plot(gUnion(polyA, polyB)).

Yet, if you select only those polygons that intersects polyB[polyA, ], dissolve will give you the expected result:

plot(gUnion(polyA, polyB[polyA, ]))

Union overlay

You can and should subset the first layer as well, gUnion(polyA[polyB, ], polyB[polyA, ]), if it has more than one feature.

edit If you want to disaggregate the multi-polygon feature into single-polygon features afterwards you can simply use the disaggregate function from the raster package.

Related Question