空间连接两个简单的特征{sf}超过1密耳.条目尽快

2023年5月18日 205次阅读

我希望这不是一件轻而易举的事情,但我真的找不到答案,而且对于这个话题来说我太新了,无法想出替代方案.所以这是问题所在：

我有两个shapefile x和y,代表Sentinel2卫星图像的不同处理级别.
x包含大约1.300.000个多边形/段,完全覆盖图像,无需任何其他重要信息.
y包含大约500个多边形,代表图像的无云区域(除了一些“云洞”外,还覆盖了大部分图像)以及有关4列中使用过的图像的信息(传感器,时间……)

即时尝试将图像信息添加到x中的位置x由y覆盖.很简单？我只是找不到一种方法来让它成为hapen而不用几天.

我将x作为一个简单的特征{sf}读取,因为使用shapefile / readOGR读取它需要很长时间.
我和y尝试过不同的事情

当我尝试合并(x,y)时,我只能拿一个sf,因为合并不支持两个sf.
合并x(作为sf)和y(作为shp)给出了错误“无法分配大小为13.0 Gb的向量”

所以我尝试了sf :: st_join(x,y),它支持两个变量为sf,但现在仍然没有完成28小时

对于10.000段子集,sf :: st_intersect(x,y)花了大约9分钟,因此对于整个片段来说可能不会快得多.

可以将x子集化为几个小块解决整个事情还是有另一个简单的解决方案？我可以对我的工作区做些什么来使合并工作,或者根本没有加入那些多边形的快捷方式？

非常感谢提前,我希望我的描述不会太模糊！

我的小工作站：

赢7比64
8 GB RAM
intel i7-4790 @ 3,6 GHz

干杯,
马蒂亚斯

最佳答案 我经常遇到这种问题,正如@ manotheshark2所承认的那样,我更喜欢在我的矢量图层中进行循环.这是我的建议：

加载您的数据

library(raster)
library(rgdal)
x <- readOGR('C:/', 'sentinelCovers')
y <- readOGR('C:/', 'cloudHoles')

指定y ID以标识哪些x多边形与y多边形相交并在x表中创建列

x$xyID <- NA # Answer col
y$yID <- 1:nrow(y@data) # ID col

运行一个子循环x

for (posX in 1:nrow(x@data)){
  pol.x <- x[posX, ]
  intX <- raster::intersect(pol.x, y)
  # x$xyID[posX] <- intX@data$yID ## Run this if there's unique y polygons
  # x$xyID[posX] <- paste0(intX@data$yID, collapse = ',') ## Run this if there's multiple y polygons
}

您可以检查是否更好地在x o y层上运行循环

x$xyID <- NA # Answer col
x$xID <- 1:nrow(x@data) # ID Col

for (posY in 1:nrow(y@data)){
  pol.y <- y[posY, ]
  intY <- tryCatch(raster::intersect(pol.y, x), finally = 'NULL')
  if (is.null(intY)) next
  x$xyID[x@data$xID %in% intY@data$xID] <- pol.y$yID
}
点赞