使用降雪和 sfLapply 在 R 中栅格化多边形答案

【问题标题】：Rasterize polygons in R using snowfall & sfLapply使用降雪和 sfLapply 在 R 中栅格化多边形
【发布时间】：2017-01-31 12:29:40
【问题描述】：

我想将一个非常大的矢量文件光栅化为 25m，并且使用 'cluster' 包取得了一些成功，调整了 qu 的 here 和 here，这对于该特定数据非常有效。

但是，我现在有一个更大的矢量文件，需要进行栅格化，并且可以访问使用降雪的集群。我不习惯集群功能，我只是不确定如何设置 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，但我就是无法让它运行。谢谢

【问题讨论】：

如果您正在寻找光栅化（大）矢量数据的速度，请查看gdalUtils::gdal_rasterize。这通常比raster::rasterize 快得多。
好的，谢谢，我也会看看
我删除了 rasFunction 并将 rDis 更改为 "rDis = sfLapply(1:n, fun=function(x) rasterize(shp[parts[[x]],], rw, 'CODE') )" 但现在我在 checkForRemoteErrors(val) 中得到错误：96 个节点产生了错误；第一个错误：“数据”必须是向量类型，为“NULL”。难住了。
@joberlin 哦。我一直在寻找一种方法来加快矢量 - > 光栅操作...
我对 gdalUtils::gdal_rasterize 的第一印象很好，本周将更新

标签： r linux cluster-computing raster rasterize

【解决方案1】：

因为我无法在 cmets 中进行格式化：

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

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

【讨论】：

不抱歉，我的意思是欢迎所有关于脚本可能出错的想法！我知道未定义和重新投影之间的区别:)

【解决方案2】：

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

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

因此，我将大型 shapefile 分段写入与处理器数量相等的部分。

最初使用 raster::rasterize 是因为我可以简单地将存储在内存中的矢量对象传递给光栅化而不会出现写入问题（尽管我希望将其写入），将此数据光栅化为 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% 用于光栅生成和光栅化）比并行 raster::rasterize 快大约 9 倍。

注意：我最初尝试将矢量分成 n 个部分，但只创建 1 个空白栅格。然后我同时从所有集群节点写入相同的空白栅格，但这破坏了栅格并使其在 R/Arc/anything 中无法使用（尽管通过函数没有错误）。上面是比较稳定的方法，但是要制作n个空白栅格而不是1个，增加了处理时间，加上合并n个栅格是额外的处理。

警告 - raster::rasterize 并行在 rasterize 函数中没有 writeRaster 而是作为单独的行，由于存储到临时文件等，这将增加原始运行中的处理时间.

编辑：为什么来自 gdal_rasterize 的栅格的频率表与 raster::rasterize 不同？我的意思是，对于 1 亿个单元，我预计会有一些差异，但对于某些代码，它有 1000 个单元不同。我以为它们都被质心光栅化了？

【讨论】：