使用 gIntersect 查找多个 SpatialPolygons 的成对交集？

Question

我想生成一个空间对象，其中包含空间多边形列表的所有成对交集。手动使用 gIntersect 处理两层的交集很容易，但我想同时找到所有成对的交集。我有 13 个多边形，所以有 156 种不同的组合来检查重叠。 lapply 或 for 循环似乎都可以工作，但我认为我需要一个包含所有可能的空间多边形组合的矩阵。

这是数据子样本的要点：

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

这是我到目前为止尝试过的方法：

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

我尝试使用 combn() 函数来修改一个有点相关的问题，其中 test 是一个空间多边形列表：

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

我收到以下错误消息：

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

以下不会产生错误，但在 ArcMap 中进行目视检查，并没有得到所有重叠的情况，可能只是最后一轮。

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

我认为我的问题是我需要 "rbind" 或将每次迭代中的所有空间多边形对象合并到一个列表中。也许如果我首先创建 SpatialPolygonDataframes 而不是 SpatialPolygons，它会更容易组合吗？我还需要避免 i 和 j 是相同的数字，否则我会得到相同的多边形自身重叠。

我尝试使用以下代码将每次迭代的重叠存储到空间对象列表中，但 spRbind 在第一个运行上阻塞，因为我认为用于存储结果的对象开始时为 NULL。

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

你的搭建方式test非常讲究。这不是一个空间对象，这使得它特别适合使用...
您使用 combn 的方式不正确。它不能使用这样的对象列表。文档说：

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

所以你应该按原样使用它。您提供的示例和数据不容易按原样使用。我不得不修改不同的东西。下次，您应该提供一个与您的数据问题相对应的简单示例。参考：How to make a great R reproducible example?
据我所了解的情况，我向您建议一种可能的方法。但是您需要了解 R 中 SpatialPolygons 对象的组织方式，以使您的脚本更清晰。

所以在这里，我以一种可用的方式使用 combn。我像您一样检索列表中的所有多边形交点。当交集不为 NULL 时，我还将其检索为 Polygons 对象，以便能够将所有多边形收集为 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")

使用 gIntersect 查找多个 SpatialPolygons 的成对交集？

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

r

spatial