3

I have a set of polygons (available here: http://ulozto.cz/x39ELbFt/polygons-zip) created in QGis and processed in R. I would like to know if they are overlapping or not, and if yes, where is the overlapped area (and between which polygons)? I know that to figure this out I can use tools:

  • over{sp} - return data.frame not a spatial object!
  • intersect{raster}
  • gIntersection{rgeos}

I can find intersection by comparing 1 by 1 polygon each time step. However as my real data are huge, I prefer to make a loop and include the iteration of polygons:

  • intersect(p1, p2)
  • intersect(p1, p3)
  • intersect(p1, p4)
  • intersect(p2, p3)
  • etc...

For the moment, I made up only loop for one polygon intersecting another 3. From resulting list of overlapping polygons I can't however really say which polygons are overlapped. And if they are not (p1 and p4) this NULL value is not included in my resulting list!

library(sp)
library(raster)
library(rgeos)
library(spatstat)
library(rgdal)     

setwd("D:/...")
p1<- readOGR(dsn=getwd(), layer="p1")
p2<- readOGR(dsn=getwd(), layer="p2")
p3<- readOGR(dsn=getwd(), layer="p3")
p4<- readOGR(dsn=getwd(), layer="p4")

ab<-list(p2, p3, p4)

# my loop only for 1 polygon, not for combinations of polygons
by two intersection methods

# intersection by gIntersection{rgeos}
int1<-list()
for (i in 1:length(ab)){
int11<-gIntersection(p1, ab[[i]], byid = FALSE)
int1[[i]]<-int11
}

# intersection by intersect{raster}
rst.int<-list()
for (i in 1:length(ab)){
rst.int1<-intersect(p1, ab[[i]])
rst.int[[i]]<-rst.int1
}

enter image description here

4

I played with assign(), ls(), and mget() to accomplish something that I believe will improve your workflow. First I use ls() to get a list of all environment variables that start with "p":

names_poly <- ls(pattern='^p.')

I used combn() to find all the unique combinations of polygons

combos <- combn(names_poly,2)

I looped over combos using get() to fetch the polygons i and j named in combos. I also used assign() to create a new object whose name is pasted together (e.g. int_p1_p2)

for(k in seq_along(combos[1,])){
    i <- combos[1,k]
    j <- combos[2,k]
    print(paste("intersecting",i,j))
    assign(paste("int",i,j,sep="_"),gIntersection(get(i), get(j), byid = FALSE))
}

You may not need this step below since you have all the new int_pi_pj SpatialPolygons as separate objects, but if you want them in a list:

int <- mget(ls(pattern='^int_.'))
rm(list = ls(pattern='^int_.')   #cleanup

Then you can access the items in your list with $:

plot(int$int_p1_p2)

EDIT: In the above, when there is no overlap, the variable will be NULL.If you don't want to even see the combinations where there is no overlap, you could add this if statement to the same for loop:

for(k in seq_along(combos[1,])){
    i <- combos[1,k]
    j <- combos[2,k]
    print(paste("intersecting",i,j))
    if(gIntersects(get(i), get(j), byid = FALSE)){
        assign(paste("int",i,j,sep="_"),gIntersection(get(i), get(j), byid = FALSE))
    }
}
  • 1
    thank you a lot! it works great and I have learn a lot about R data processing and its almost immense possibilities ! Great idea to use combn() and another tools ! only one remark: rm(ls(pattern='^int_.') #cleanup doesn't work for my, but rm(list = ls(pattern='^int_.')) does. – maycca Aug 3 '15 at 2:44
  • My pleasure! Thinking about this also helped me with my own work so everyone wins! I edited the answer to include the list = in the rm() function. I'm not sure why some systems need the argument list specific, but it does work on mine so including it makes it a more general answer – Remi Daigle Aug 3 '15 at 16:56

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.