Loop to check multiple polygons overlap in r

Question

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
}

Answer 1

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))
    }
}