폭설을 사용하여 R의 폴리곤을 래스터 화합니다.

Question

매우 큰 벡터 파일을 25m로 래스터 화하고 '클러스터'패키지로 성공한 적이 있습니다. 특정 데이터에 대해 잘 작동하는 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 

내 전체 코드 : 나는 기능과 sfLapply을 변경, 여러 가지를 시도했습니다

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() 

, 하지만 난 그냥 달릴 수 없어. 감사

Answer 1

은 그래서 눈을 포기하고 나는 gdalUtils :: gdal_rasterize에 대신 고개를 사용하여 혜택을 많이 발견

상황 & 문제을 (사람이 대답 할 수있을 것으로 하나의 단점으로?) 내 벡터 데이터는 ESRI 파일 Geodatabase 내부에 존재하며 일부 사전 사전 래스터 화가 필요합니다. 문제는 없지만 rgdal :: readOGR이 좋습니다. 그러나 gdal_rasterize는 벡터 데이터에 대한 경로 이름이 필요하기 때문에 처리 된 벡터 데이터를 쓸 수 없기 때문에 문제가 발생했습니다. 지오 데이터베이스 외부의 셰이프 파일의 최대 파일 크기를 초과하며 gdal_rasterize는 개체를 허용하지 않으며 .gdbs 또는 .Rdata/.rds 파일. 개체를 gdal_rasterize에 전달하는 방법 ??

그래서 대형 셰이프 파일을 프로세서 수와 동일한 세그먼트로 작성했습니다.

원래 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) 

시간이 코드의 전체 작업 (래스터 생성 및 래스터 화용 & 40 %를 기록/벡터 판독/처리하고, 60 %의 분리)의 주위에 평행 한 래스터 :: 래스터 화보다 더 빠른 9 번이었다.

참고 : 처음에는 벡터를 n 부분으로 분할했지만 빈 래스터는 1 개만 작성하여이 작업을 시도했습니다. 그런 다음 모든 클러스터 노드에서 동일한 빈 래스터에 동시에 썼습니다.하지만 래스터가 손상되어 R/Arc/anything에서 오류없이 사용할 수 없게되었습니다. 위의 방법은보다 안정적이지만, 1 대신 빈 n 개의 래스터를 만들어야하므로 처리 시간이 늘어나고 n 개의 래스터를 병합하는 것은 별도의 처리입니다.

주의 - 래스터 :: 병렬 래스터는 래스터의 함수 내에서가 아니라 임시 파일 등

로 인해 스토리지에 원래 실행에서 처리 시간이 증가하는 별도의 행으로 writeRaster이 없었다

EDIT : gdal_rasterize의 래스터의 빈도 테이블이 raster :: rasterize와 다른 이유는 무엇입니까? 내 말은 1 억 개의 셀이 있는데 약간의 차이는 있지만 몇 가지 코드의 경우 1000 개의 셀이 다른 것을 의미합니다. 둘 다 중심에서 래스터 화 한 줄 알았는데?

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 

가 투사 덮어 쓰기 = 다시 투영 데이터입니다. 좋아