两个简单特征 {sf} 上的空间连接超过 1 百万。尽快输入

Question

我希望这不是太琐碎，但我真的找不到答案，而且我对这个话题太陌生，无法自己提出替代方案。所以这里是问题：

我有两个 shapefile x 和 y，它们代表 Sentinel2 卫星图像的不同处理级别。x 包含大约 1.300.000 个多边形/段，完全覆盖了图像扩展，没有任何进一步的重要信息。
y 包含大约 500 个多边形，代表图像的无云区域（也覆盖了图像的大部分，除了一些“云洞”）以及 4 列（传感器、时间...）

我正在尝试将图像信息添加到 x 中 x 被 y 覆盖的位置。很简单？我就是找不到不花几天时间就可以实现的方法。

我将 x 作为一个简单的特征 {sf} 读入，因为使用 shapefile / readOGR 读取它需要很长时间。我用 y 尝试了不同的东西

当我尝试合并（x，y）时，我只能取一个 sf，因为合并不支持两个 sf。合并 x（作为 sf）和 y（作为 shp）给我错误“无法分配大小为 13.0 Gb 的向量”

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

sf::st_intersect(x,y) 10.000 段子集大约需要 9 分钟，因此对于整个片段来说可能不会快很多。

可以将 x 子集为几个较小的部分来解决整个问题还是有另一个简单的解决方案？我可以对我的工作区做些什么来使合并工作，还是根本没有加入那么多多边形的捷径？

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

我的小工作站：
win 7 64 位
8 GB RAM
intel i7-4790 @ 3,6 GHz

Answer 1

我经常遇到这类问题，正如 @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
}

您可以检查在 xoy 层上运行循环是否更好

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
}