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

2023年5月10日 313次阅读

我有一个包含15,000个空间线的巨大数据集,我使用37,000个点的所有组合创建了这些数据集.对于每个空间线,我想提取线接触的多边形(或栅格 – 更快)的最大值.从本质上讲,这是Arc术语中非常大的“空间连接”.如果在多边形图层上覆盖线条,则输出将是所有属性字段中空间线条的最大值 – 每个字段代表一年中的一个月.我已经包含了一个栅格数据集,该栅格数据集仅在1990年1月创建的多边形文件的分辨率为~30m – 栅格代表了一种我认为可以节省时间的替代方法.多边形&栅格图层代表一个大的空间区域：大约30km×10km.数据可用
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的连续变量运行此代码.

最佳答案 这是一个应该给出一个很好的近似的黑客.它可能会得到改善(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 /长度(行))/(3600 * 24)约50天……

50台电脑只用了1天

请注意,使用更多行可以获得相对更高的效率(因为要查询的唯一单元相对较少).

点赞