5

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.

  • Some sample data with both overlaying and overlaying polygons would be helpfull. – dof1985 Aug 13 '16 at 10:45
1

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.

  • Thanks for your comprehensive answer. I will try to see if I can implement it on my data and write back. – joaoal Aug 14 '16 at 7:30
  • @joaoal if it works you can mark it as answer to make it more available to other members of the forum – dof1985 Aug 14 '16 at 12:27
  • i put the solution above. Thanks for your help. I'm unfortunately not allowed yet to upvote. – joaoal Aug 15 '16 at 10:28
  • 1
    If I am not mistaken, your code creates a single feature (multipolygon) out of all overlapping ones. In your example it does not appear because only two polygons are unioned together. However I would prefer multiple features each containing the overlapping polygons. A solution could be to dissagregate afterwards like indicated here: stackoverflow.com/questions/19315832/… – joaoal Aug 15 '16 at 10:43
  • 1
    it worked. I edited your answer above slightly adding this detail. – joaoal Aug 15 '16 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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