使用 snowfall 和 sfLapply 栅格化 R 中的多边形

Question

我想将一个非常大的矢量文件光栅化到25m，并且在"cluster"包方面取得了一些成功，在这里和这里调整了qu，这对特定的数据非常有效。

然而，我现在有一个更大的矢量文件，需要光栅化，并且可以访问使用降雪的集群。我不习惯集群函数，只是不知道如何设置sfLapply。在集群中调用sfLapply时，我一直收到以下错误：

Error in checkForRemoteErrors(val) : 
one node produced an error: 'quote(96)' is not a function, character or symbol
Calls: sfLapply ... clusterApply -> staticClusterApply -> checkForRemoteErrors

我的完整代码：

library(snowfall)
library(rgeos)
library(maptools)
library(raster)
library(sp)
setwd("/home/dir/")
# Initialise the cluster...
hosts = as.character(read.table(Sys.getenv('PBS_NODEFILE'),header=FALSE)[,1]) # read the nodes to use
sfSetMaxCPUs(length(hosts)) # make sure the maximum allowed number of CPUs matches the number of hosts
sfInit(parallel=TRUE, type="SOCK", socketHosts=hosts, cpus=length(hosts), useRscript=TRUE) # initialise a socket cluster session with the named nodes
sfLibrary(snowfall)
# read in required data
shp <- readShapePoly("my_data.shp")
BNG <- "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"
crs(shp) <- BNG
### rasterize the uniques to 25m and write (GB and clipped) ###
rw <- raster(res=c(25,25), xmn=0, xmx=600000, ymn=0, ymx=1000000, crs=BNG)
# Number of polygons features in SPDF
features <- 1:nrow(shp[,])
# Split features in n parts
n <- 96
parts <- split(features, cut(features, n))
rasFunction = function(X, shape, raster, nparts){
ras = rasterize(shape[nparts[[X]],], raster, 'CODE')
return(ras)
}
# Export everything in the workspace onto the cluster...
sfExportAll()
# Distribute calculation across the cluster nodes...
rDis = sfLapply(n, fun=rasFunction,X=n, shape=shp, raster=rw, nparts=parts) # equivalent of sapply
rMerge <- do.call(merge, rDis)
writeRaster(rMerge, filename="my_data_25m",  format="GTiff", overwrite=TRUE)
# Stop the cluster...
sfStop()

我尝试了很多方法，更改函数和sfLapply，但我就是无法运行。感谢

Answer 1

因为我不能在注释中进行格式化：

library(maptools)
shp <- readShapePoly("my_data.shp")
BNG <- "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"
shp.2 <- spTransform(shp, BNG)
#Continue as before

覆盖投影！=重新投影数据。

Answer 2

好吧，所以我放弃了降雪，转而研究了gdalUtils:gdal_rasterize，发现使用它有很多好处(有一个缺点，有人可能会回答？)

上下文&问题：我的矢量数据存在于ESRI文件Geodatabase中，需要进行一些预光栅化处理。没问题，rgdal:：readOGR很好。然而，由于gdal_rasterize需要矢量数据的路径名，我在这里遇到了麻烦，因为我无法写出处理过的矢量数据，它们超过了地理数据库外shapefile的最大文件大小，gdal_resterize将不接受对象、.gdb或.Rdata/.rds文件的路径如何将对象传递给gdal_rasterize

所以我把这个大的shapefile写在与处理器数量相等的段中。

最初使用光栅化：：光栅化是因为我可以简单地将存储在内存中的矢量对象传递给光栅化，而不会出现写入问题(尽管我希望写入它)，将这些数据光栅化到25m。这花了相当长的时间，即使是并行的。

解决方案：并行的gdal_rasterize。

# gdal_rasterize in parallel
require(gdalUtils)
require(rgdal)
require(rgeos)
require(cluster)
require(parallel)
require(raster)
# read in vector data
shape <- readOGR("./mygdb.gdb", layer="mydata",stringsAsFactors=F)
## do all the vector processing etc ##
# split vector data into n parts, the same as number of processors (minus 1)
npar <- detectCores() - 1
features <- 1:nrow(shape[,])
parts <- split(features, cut(features, npar))
# write the vector parts out
for(n in 1:npar){
writeOGR(shape[parts[[n]],], ".\parts", paste0("mydata_p",n), driver="ESRI Shapefile")
}
# set up and write a blank raster for gdal_rasterize for EACH vector segment created above
r <- raster(res=c(25,25), xmn=234000, xmx=261000, ymn=229000, ymx=256000, crs=projection(shape))    
for(n in 1:npar){
writeRaster(r, filename=paste0(".\gdal_p",n,".tif"), format="GTiff", overwrite=TRUE)
}
# set up cluster and pass required packages and objects to cluster
cl <- makeCluster(npar)
clusterEvalQ(cl, sapply(c('raster', 'gdalUtils',"rgdal"), require, char=TRUE))
clusterExport(cl, list("r","npar"))
# parallel apply the gdal_rasterize function against the vector parts that were written, 
# same number as processors, against the pre-prepared rasters
parLapply(cl = cl, X = 1:npar, fun = function(x) gdal_rasterize(src_datasource=paste0(".\parts\mydata_p",x,".shp"),
dst_filename=paste0(".\gdal_p",n,".tif"),b=1,a="code",verbose=F,output_Raster=T))
# There are now n rasters representing the n segments of the original vector file
# read in the rasters as a list, merge and write to a new tif. 
s <- lapply(X=1:npar, function(x) raster(paste0(".\gdal_p",n,".tif")))
s$filename <- "myras_final.tif"
do.call(merge,s)
stopCluster(cl)

此代码中整个作业的时间(矢量读取/处理/写入的时间为60%，光栅生成和光栅化的时间为40%)大约是并行光栅化的9倍。

注意：我最初尝试将矢量拆分为n个部分，但只创建了一个空白光栅。然后，我从所有集群节点同时写入同一个空白光栅，但这损坏了光栅，使其在R/Arc/任何东西中都不可用(尽管在执行函数时没有出错)。以上是一种更稳定的方法，但必须制作n个空白光栅而不是1个，这增加了处理时间，再加上合并n个光栅是额外的处理。

警告-光栅：：并行光栅化在光栅化函数中没有writeRaster，而是作为一个单独的行，由于存储到临时文件等，这将增加原始运行中的处理时间。

EDIT：为什么gdal_rasterize中光栅的频率表与rasterize:：rasterize不同？我的意思是，对于1亿个单元格，我预计会有一些不同，但对于某些代码，它有几千个单元格不同。我以为它们都是通过质心光栅化的？