在R中提取栅格的最快方法(改善我的可重复代码的时间)
在R中提取栅格的最快方法(改善我的可重复代码的时间)
提问于 2018-06-15 14:18:08
我想知道我是否已经最大化了可以提取栅格点周围缓冲区域的平均值的速度。
还能在本地进一步提高性能吗？我已经在使用parallel mclapply了，我知道我可以通过在集群上设置和运行它来获得进一步的收益(使用集群或获得更多的cpu不是我想要的答案)。
复制一些数据：
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 rastertest2：，但是如果必须从VR生成velox (模拟具有不同栅格)，则具有velox rastertest3：，但是使buffertest4：必须生成两者(从VR)test5：生成velox，从tif生成velox，如果具有buffertest6：，则生成两者(tif版本)
#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并行运行，您的结果会因内核而异)：
#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有没有人能给我提点建议，让它变得更快，或者我的本地速度已经达到极限了？，谢谢！
回答 1
Stack Overflow用户
发布于 2018-06-16 07:02:02
我从@dbaston得到了一个建议，在velox中预先裁剪栅格。这是到目前为止我发现的在R中提取栅格的最快方法：
如果你已经有了velox栅格，这是难以置信的快(闪电)，即使你必须在函数中加载缓冲区(没有显示，但在我的系统上大约0.4个小时过去了)：
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 即使您必须动态加载不同的栅格，速度也很快(多次模拟加载相同的栅格)：
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 页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/50870080
复制相关文章
相似问题