Растеризация полигонов в R с использованием снегопада и sfLapply

Question

Я хотел бы растеризовать очень большой векторный файл до 25 м и добился некоторого успеха с пакетом '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

Итак, я отказался от снегопада и вместо этого заглянул в gdalUtils::gdal_rasterize и обнаружил множество преимуществ его использования (с одним недостатком, на который кто-то мог бы ответить?)

Контекст и проблема: Мои векторные данные существуют в файловой базе геоданных ESRI и требуют некоторой предварительной обработки растеризации. Нет проблем, rgdal::readOGR в порядке. Однако, поскольку gdal_rasterize требует путь к векторным данным, у меня возникли проблемы, поскольку я не мог записать обработанные векторные данные, они превышают максимальный размер файла для шейп-файла вне базы геоданных, и gdal_rasterize не будет принимать объекты, пути к.gdbs или.Rdata/.rds файлы. Как передать объект в gdal_rasterize??

Поэтому я выписал большой шейп-файл в сегментах, равных количеству процессоров.

Первоначально использовался raster::rasterize, поскольку я мог просто передать растровый векторный объект, сохраненный в памяти, для растеризации без проблем с записью (хотя мне бы хотелось, чтобы он был записан), растеризовав эти данные до 25 м. Это заняло довольно много времени, даже параллельно.

Решение: 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 частей, но создав только 1 пустой растр. Затем я одновременно записал в один и тот же пустой растр со всех узлов кластера, но это повредило растр и сделало его непригодным для использования в R/Arc/ что угодно (несмотря на то, что функция прошла без ошибок). Выше приведен более стабильный способ, но вместо 1 необходимо создать n пустых растров, что увеличивает время обработки, плюс объединение n растров - дополнительная обработка.

У caveat - raster::rasterize параллельно не было writeRaster внутри функции rasterize, а скорее как отдельная строка, которая будет увеличивать время обработки при первоначальном запуске из-за хранения временных файлов и т. д.

РЕДАКТИРОВАТЬ: Почему таблицы частот из растра из gdal_rasterize не совпадают с растром::rasterize? Я имею в виду, что при 100 миллионах ячеек я ожидаю небольшую разницу, но для некоторых кодов это было несколько разных ячеек. Я думал, что они оба растеризованы по центроиду?

Answer 2

Потому что я не могу сделать форматирование в комментариях:

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

Перезапись проекции!= Перепроектирование данных.