在R中提取栅格的最快方法(缩短了可复制代码的时间) [英] Fastest way to extract a raster in R (improve the time of my reproducible code)

问题描述

我想知道我是否已最大化提取栅格点周围缓冲区域的均值的速度.

I'm wondering if I have maximized the speed at which a mean of an area buffered around a point in a raster can be extracted.

可以在这些位置上进一步提高性能吗?我已经使用了parallel mclapply，而且我知道我可以通过在群集上进行设置和运行来获得进一步的收益(使用群集或获得更多的cpu并不是我要找的答案).

Can performance be improved any further on these LOCALLY? I use parallel mclapply already, and I know I could get further gains by setting up and running this on a cluster (use a cluster or get more cpu's is not the answer I'm looking for).

复制一些数据:

library(raster)
library(parallel)
library(truncnorm)
library(gdalUtils)
library(velox)
library(sf)
ras <- raster(ncol=1000, nrow=1000, xmn=2001476, xmx=11519096, ymn=9087279, ymx=17080719)
ras[]=rtruncnorm(n=ncell(ras),a=0, b=10, mean=5, sd=2)
crs(ras) <- "+proj=utm +zone=51 ellps=WGS84"

writeRaster(ras,"testras_so.tif", format="GTiff")

gdalbuildvrt(gdalfile = "testras_so.tif", 
             output.vrt = "testvrt_so.vrt")

x1 <- runif(100,2001476,11519096)
y1 <- runif(100, 9087279,17080719)

poly <- st_buffer(st_sfc(st_point(c(x1[1],y1[1]), dim="XY"),crs=32651),200000)
vras <- velox("testvrt_so.vrt")
###########

测试:

• test1:如果具有多边形和velox栅格
• test2:如果必须生成缓冲区但具有velox栅格
• test3:如果必须从VR中生成velox(模拟具有不同的栅格)但具有缓冲区
• test4:必须同时生成两者(通过VR)
• test5:从具有缓冲区的tif生成velox
• test6:同时生成(tif版本)

• test1: if have poly and velox raster
• test2: if have to generate buffer but have velox raster
• test3: if have to generate velox from VR (simulating having different rasters) but having the buffer
• test4: have to generate both (from VR)
• test5: generate velox from tif have buffer
• test6: generate both (tif version)

#Test time if have poly and velox raster
test1 <- system.time(mclapply(seq_along(x1), function(x){
  vras$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate buffer but have velox raster
test2 <- system.time(mclapply(seq_along(x1), function(x){
  vras$extract(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000), fun=function(t) mean(t,na.rm=T))
}))


#Test time if have to generate velox from VR (simulating having different rasters) but having the buffer
test3 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testvrt_so.vrt")$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from VR AND generate buffer (simulating a list of rasters with different buffers each)
test4 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testvrt_so.vrt")$extract(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000), fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from TIF (simulating having different rasters) but having the buffer
test5 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testras_so.tif")$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from TIF AND generate buffer (simulating a list of rasters with different buffers each)
test6 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testras_so.tif")$extract(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000), fun=function(t) mean(t,na.rm=T))
}))

我的结果(由于mclapply并行运行，您的内核会有所不同):

My results (yours will vary with cores due to mclapply running parallel):

#Test time if have poly and velox raster
   > test1
   user  system elapsed 
  0.007   0.022   3.417 

#Test time if have to generate buffer but have velox raster
> test2
   user  system elapsed 
  0.007   0.023   3.540 

#Test time if have to generate velox from VR (simulating having different rasters) but having the buffer
> test3
   user  system elapsed 
  0.016   0.037  10.366 

#Test time if have to generate velox from VR AND generate buffer (simulating a list of rasters with different buffers each)
> test4
   user  system elapsed 
  0.017   0.035  10.309 

#Test time if have to generate velox from TIF (simulating having different rasters) but having the buffer
> test5
   user  system elapsed 
  0.015   0.033   9.258 

#Test time if have to generate velox from TIF AND generate buffer (simulating a list of rasters with different buffers each)
> test6
   user  system elapsed 
  0.016   0.034   9.347

有人可以提出任何建议来加快速度吗?还是我在这里达到了本地速度?谢谢！

推荐答案

我从

I got a suggestion to pre-crop the raster in velox from @dbaston. This is the fastest way I have found to extract raster so far in R:

如果您已经有了velox栅格，即使您必须将缓冲区加载到函数中(未显示，但在我的系统上大约0.4时它就已经出来了)，这都非常快(闪电):

If you have the velox raster already this is unbelievably fast (lightning), even if you have to load the buffer in the function (not shown, but it comes out on my system around 0.4 elapsed):

test7_lightning <- system.time(mclapply(seq_along(x1), function(x){
  q <- vras$crop(poly);vras$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

> test7_lightning
   user  system elapsed 
  0.001   0.005   0.355 

即使您必须动态加载不同的栅格(模拟多次加载相同的栅格)也可以快速完成:

Fast even if you have to dynamically load different rasters (simulated loading same raster many times):

test8 <- system.time(mclapply(seq_along(x1), function(x){ ras<-velox("testras_so.tif");ras$crop(poly);ras$extract(poly, fun=function(t) mean(t,na.rm=T)) }))

test9 <- system.time(mclapply(seq_along(x1), function(x){ ras<-velox("testras_so.tif");ras$crop(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000));ras$extract(poly, fun=function(t) mean(t,na.rm=T)) }))


> test8
   user  system elapsed 
  0.011   0.016   4.450 
> test9
   user  system elapsed 
  0.006   0.012   4.333 

这篇关于在R中提取栅格的最快方法(缩短了可复制代码的时间)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！