从空间线更有效地覆盖多边形或 extract() 栅格数据答案

【问题标题】：More efficiently overlay polygon or extract() raster data from spatial lines从空间线更有效地覆盖多边形或 extract() 栅格数据
【发布时间】：2015-09-16 11:08:44
【问题描述】：

我使用 37,000 个点的所有组合创建了一个包含 15 亿条空间线的庞大数据集。 对于每条空间线，我想提取该线所接触的多边形（或栅格 - 任何更快的）的最大值。本质上这是 Arc 术语中的一个非常大的“空间连接”。如果在多边形图层上叠加线，则输出将是所有属性字段中空间线的最大值 - 每个字段代表一年中的一个月。我还包含了一个栅格数据集，它是从 1990 年 1 月创建的多边形文件，分辨率约为 30m - 栅格代表了一种我认为可以节省时间的替代方法。多边形和栅格图层代表了一个很大的空间区域：大约 30 公里 x 10 公里。数据可用here。我在 .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”似乎并没有节省时间。 还有其他方法可以更有效地完成这种可怕的叠加或提取（）吗？请注意，这些数据目前被分类为“1”或“0”，但最终我想运行这段代码以 0:300 运行的连续变量。

【问题讨论】：

所有 37,000 点对（不包括零长度线 A-A）应该只给您 684,481,500 行来检查，因为 A-B 与 B-A 命中相同的多边形。所以这是 422 天...

标签： r overlay extract gis

【解决方案1】：

这是一个hack，应该给出一个很好的近似值。它可能会得到改进（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 天

请注意，越多行越高效（因为要查询的唯一单元格相对较少）。

【讨论】：

【解决方案2】：

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

但是我不会在 R 中对它们进行栅格化。我会编写一些 C 代码来获取栅格和 37,000 个点位置的数据，然后使用 Bresenham 画线算法来获取线的栅格位置.在这些位置对栅格进行采样，然后对这些数据执行任何您需要的操作。 Bresenham 算法的快速代码应该很容易获得，你甚至可以找到在 GPU 上运行的版本以实现巨大的加速。什么东西画直线比显卡快？

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

或者从亚马逊（或其他云提供商）租用 1000 台服务器半天。

【讨论】：