更高效地叠加多边形或从空间线提取（extract()）栅格数据

Question

我有一个巨大的数据集，由37000个点的所有组合生成了15亿个空间线。对于每个空间线，我想提取它所接触到的多边形（或栅格 - 以更快的方式为准）的最大值。本质上，这是一个非常大的“空间连接”在Arc术语中。如果将线覆盖在多边形图层上，则输出将是所有属性字段的空间线的最大值 - 每个字段代表一年中一个月。我还包括了一个栅格数据集，它是从多边形文件1990年1月份的数据以约30m分辨率创建的 - 栅格代表了我认为可以节省时间的另一种方法。多边形和栅格图层表示一个大的空间区域：大约30km x 10km。数据可在此处获得。我在.zip中包含的空间线数据集仅有9900条线，是从15亿条线的整个数据集中随机抽样的。

首先读入数据。

#polygons

 poly<-readShapePoly("ls_polys_bin",proj4string=CRS("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"))
 poly$SP_ID<-NULL #deleting this extra field in prep for overlay

#raster - this represents only one month (january 1990)
   #raster created from polygon layer but one month only

     raster.jan90<-readGDAL("rast_jan90.tif") 
     raster.jan90<-raster(raster.jan90) #makes it into a raster

#lines (9900 of 1.5 billion included)

     lines<-readShapeLines("l_spatial",proj4string=CRS("+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"))

为了更方便管理行数据，取样50行。

 lines.50<-lines[sample(nrow(lines),50),]

将三个图层一起绘制

plot(raster.jan90)#where green=1
plot(poly, axes=T,cex.axis=0.75, add=T)
plot(lines.50, col="red", add=TRUE)

首先，我尝试使用叠加，但按照当前速度，在我的机器上运行全部数据集的15亿条数据需要大约844天。

 ptm <- proc.time() #start clock
 overlays.all<-over(lines.50,poly, fn=max)
 ptm.sec.overlay<-proc.time() - ptm # stop clock
 ptm.sec.overlay #.56 sec w/ n=12 lines; 2.3 sec w/ 50 lines

接下来，我将多边形转换为栅格图像（仅限1990年1月），并使用空间线运行了一个extract()，但这需要更长的时间。

 ptm <- proc.time() # Start clock
 ext.rast.jan90<-extract(raster.jan90,lines.50, fun=max, method=simple)
 ptm.sec.ext<-proc.time() - ptm # stop clock
 ptm.sec.ext #32 sec w/ n=12 lines; 191 sec w/ n=50 lines

我的尝试将所有“0”单元格转换为“NA”似乎没有节省时间。是否有另一种更高效的方法来处理这个庞大的覆盖或提取（extract()）？请注意，这些数据目前被分为“1”或“0”，但最终我想对一个连续变量运行0:300的代码。

- Luke

---

所有37,000个点对（不计零长度的A-A线）应该只会生成684,481,500条待检查的直线，因为A-B和B-A会撞在同一个多边形上。所以大概需要422天…… - Spacedman

Answer 1

我认为最快的方法是将线条光栅化到与栅格数据相同的光栅中。

但是我不会在R中进行光栅化。我会编写一些C代码，该代码获取栅格和37,000个点位置的数据，然后使用Bresenham线绘制算法获取线的栅格位置。在这些位置对栅格进行采样，并对该数据执行所需操作。快速的Bresenham算法代码应该很容易获取，甚至可以找到用于GPU的版本以获得大规模加速。有什么比图形卡更快地绘制直线的呢？

我假设您的空间线是两个点之间的单条直线段。

或者只需租用亚马逊（或其他云提供商）的1000台服务器半天即可。

- Spacedman

Answer 2

这里有一个技巧，可以给出一个很好的近似值。它可能可以改进（getCrds需要很长时间），包括采取更大的步骤（无论您是否同意，我不知道）。

library(raster)
raster.jan90 <- raster("rast_jan90.tif") 
lines <- shapefile("l_spatial.shp", p4s="+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs")  
lines.50<-lines[sample(nrow(lines),50),]

test <- function(lns) {

  getCrds <- function(i) {
    p <- z[[i]][[1]]
    s <- (p[2,] - p[1,]) / res(raster.jan90)
    step <- round(max(abs(s)))
    if ( step < 1 ) {
        # these probably should not exist, but they do
        return( cbind(i, cellFromXY(raster.jan90, p[1, , drop=FALSE])) )
    }
    x <- seq(p[1,1], p[2,1], length.out=step)
    y <- seq(p[1,2], p[2,2], length.out=step)
    cbind(i, unique(cellFromXY(raster.jan90, cbind(x, y))))
  }

  z <- coordinates(lns)
  crd <- sapply(1:length(z), getCrds )
  crd <- do.call(rbind, crd)

  e <- extract(raster.jan90, crd[, 2])
  tapply(e, crd[,1], max)
}

system.time(res <- test(lines.50))
#  user  system elapsed 
#  0.53    0.01    0.55 

system.time(res <- test(lines))
#  user  system elapsed 
#  59.72    0.85   60.58 

(684481500 * 60.58 / length(lines)) / (3600 * 24) 大约为50天...

在50台计算机上只需要1天

请注意，随着行数的增加，效率会相对更高（因为要查询的唯一单元格相对较少）。

- Robert Hijmans