1

我想将一个非常大的矢量文件光栅化为 25m，并在 'cluster' 包中取得了一些成功，调整了 qu's 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，但我就是无法让它运行。谢谢

4

2 回答 2

0

因为我不能在评论中进行格式化：

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

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

于 2017-02-02T20:54:32.270 回答
0

好的，所以我放弃了降雪，而是研究了 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 个单元不同。我以为他们都被质心光栅化了？

于 2017-02-03T14:40:31.980 回答