using gIntersect to find pair-wise intersection of multiple SpatialPolygons?

Question

I'd like to produce a spatial object containing all the pairwise intersections for a list of spatial polygons. It's easy enough to manually use gIntersect for the intersection of just two layers, but I'd like to find all pair-wise intersections simultaneously. I have 13 polygons so there are 156 different combinations of pairs to check for overlap. Either lapply or a for loop seems like it might be able to work, but I think I'd need a matrix of all the possible combinations of spatial polygons.

Here is a gist for a subsample of the data:

https://gist.github.com/dwwolfson/b1dc7b9c084233a4a36401f7e7061897

Here is what I've tried so far:

# Overlap in spring and summer of 2016
library(dplyr)
library(gdata)
library(sp)
library(rgdal)
library(rgeos)
library(raster)
library(spatstat)

#Import Data  (from gist)
df16<-df16.sub
#Extract only locations between April 1 to August 1, 2016
sub16<-df16[df16$loctime>"2016-03-31"&df16$loctime<"2016-08-01",]
#Subset by age
colts.16<-sub16[sub16$age=="colt",]
df<-colts.16
df$id<-as.factor(df$id)

#extract coordinates
coordinates(df)<-c("location.long", "location.lat")
#define coordinate ref system
proj4string(df)<-CRS("+proj=longlat +datum=WGS84")
#reproject
df<-spTransform(df, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5
                        +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "))
colnames(df@coords)<-c("x", "y")
sp.df<-as.data.frame(df)
crane.names<-as.list(unique(sp.df$id))
testdat<-sp.df

temp.crane<-NULL
temp.crane.day<-NULL
dat.summary<-NULL
tempdat<-NULL
firstday<-NULL
lastday<-NULL
dat.summary<-NULL
i<-NULL
j<-NULL
#loop to condense points into just roost sites
for(i in 1:length(unique(testdat$id))){
  temp.crane<-testdat[testdat$id==crane.names[[i]],]
  if(nrow(temp.crane)>0){
    current.crane<-as.character(crane.names[[i]])
    temp.days<-as.list(unique(temp.crane$day))
    days.vec<-unlist(temp.days)
    firstday<-days.vec[1]
    lastday<-last(days.vec)
    for(j in seq.int(from=firstday, to=lastday)){
      tempdat<-temp.crane[which(temp.crane$day==j),]    
      firstrow<-tempdat[1,]
      lastrow<-tempdat[nrow(tempdat),]
      daylocs<-rbind(firstrow, lastrow)

      if((j-1)<firstday)
      {daylocs$night.dist<-NA}
      else{
        prevdat<-temp.crane[which(temp.crane$day==j-1),]    
        firstrow<-prevdat[1,]
        lastrow<-prevdat[nrow(prevdat),]
        prevdaylocs<-rbind(firstrow, lastrow)

        daylocs$night.dist<-sqrt((daylocs$x[1]-prevdaylocs$x[2])^2+(daylocs$y[1]-prevdaylocs$y[2])^2)
      }
            dat.summary<-rbind(dat.summary, daylocs)
          }
  }
}

roosts<-dat.summary
temp.buffer<-NULL
temp.crane<-NULL
current.crane<-NULL
temp.id<-NULL
temp.sp.df<-NULL
buff.spdf<-NULL
crane.names<-as.list(unique(roosts$id))
test<-NULL
#create list of SpatialPolygons
for(i in 1:length(unique(roosts$id))){
  temp.crane<-roosts[roosts$id==crane.names[[i]],]
  current.crane<-as.character(crane.names[[i]])
  coordinates(temp.crane)<-c("x","y")
  proj4string(temp.crane)<-CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5
                          +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")
  temp.buffer<-gBuffer(temp.crane, width=3000)
  test[i]<-list(temp.buffer)
}

I tried to adapt a question about something that kinda related using the combn() function where test is a list of spatial polygons:

combos1<-combn(test,2)
for(k in length(combos1)){
i<-combos[1,k]
j<-combos[2,k]
print(paste("intersecting", i, j))
int[k]<-gIntersection(i,j, byid=T)
}

I get the following error message:

Error in identical(spgeom1@proj4string, spgeom2@proj4string) :                         
trying to get slot "proj4string" from an object of a basic class ("list") with no slots

The following doesn't produce an error, but on visual inspection in ArcMap, it's not getting all the occasions of overlap, probably just the last round.

   overlap<-list()
for(i in length(test)){
for(j in length(test)){
overlap<-gIntersection(test[[i]],test[[j]])
} 
}

I think my problem is that I need to "rbind" or merge all the spatial polygon objects from each iteration into one list. Maybe if I first create SpatialPolygonDataframes instead of SpatialPolygons it'd be easier to combine? I also need to avoid i and j being the same number or else I'll get the same polygon overlapping itself.

I tried to store each iteration's overlap into a list of spatial objects using the following code, but the spRbind chokes on the first run because object intended to store the results starts out as NULL, I think.

library(maptools)
over.list<-NULL

for(i in length(test)){
for(j in length(test)){
overlap<-gIntersection(test[[i]],test[[j]])
over.list<-spRbind(over.list, overlap)
} 
}

Answer 1

The way you build test is very particular. This is not a Spatial object, which makes it particular to use...
The way you use combn is not correct. It can not use a list of objects like this. The documentation says:

?combn
x vector source for combinations, or integer n for x <- seq_len(n).

So you should use it as is. The example and data you give is not easy to be used as is. I had to modify differents things. Next time, you should provide a simple example that corresponds to your data problem. Refer to: How to make a great R reproducible example?
As far as I can understand from what your looking at, I propose you a possible way. But you need to learn what is the organization of SpatialPolygons objects in R to make your script cleaner.

So here, I use combn in a usable way. I retrieve all polygons intersections in a list as you did. When the intersection is not NULL, I also retrieve it as a Polygons object to be able to gather all Polygons as a SpatialPolygons:

combos <- combn(length(test),2)
int <- vector("list", length = ncol(combos))
Polygons.list <- list()
k.poly <- 0 # Number of non-null polygons
for(k in 1:ncol(combos)){
  i<-combos[1,k]
  j<-combos[2,k]
  print(paste("intersecting", i, j))
  int.tmp <- gIntersection(test[[i]],test[[j]], byid=FALSE)
  if (!is.null(int.tmp)) {
    k.poly <- k.poly + 1
    int[[k]] <- int.tmp
    Polygons.list[[k.poly]] <- int[[k]]@polygons[[1]]
    Polygons.list[[k.poly]]@ID <- paste(k.poly)
  }
}

all.Intersections <- SpatialPolygons(Polygons.list)
plot(all.Intersections, col = "red")