6

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

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

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

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

当我尝试合并（x，y）时，我只能取一个 sf，因为合并不支持两个 sf。合并 x (as sf) 和 y (as 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

4

1 回答 1

2

我经常遇到这类问题，正如@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
}
于 2017-04-05T16:10:12.017 回答