在R中将rasterstack转换为csv并行处理

Question

我有大型栅格 (>20 GB)，并希望将每个栅格转换为特殊格式的 csv 文件，如下所示：

unique_key_column
x_coordinate
y_coordinate
layer1_values
layer2_values

等等

library(raster)
r <- raster(nrows=10,ncols=10)
r[] <- rnorm(10)
stack <- stack(r,r,r,r,r)

   #create function to convert coordinate to special format
   # -34.9 will be 1034900000 
   # sxxxdddddd, where s= sign (-=1, +=2), x=degrees (34=034), 
   # and d = decimal (.9=900000)

formatCoordinate <- function(x){
  first_part <- ifelse(x < 0 , "1","2")
  second_part <- abs(as.integer(x))
  #make sure 3 part has 6 decimal places, then convert it to string
  third_part <- substr(gsub(".+\\.","",as.character(format(round(x, 2),
                             nsmall = 6))),1,6)
  result <- sprintf("%s%03d%s",first_part,second_part,third_part)
  result
}

  #the actual processing

stack =readAll(stack)
names(stack) <-c("l1", "l2", "l3", "l4", "l5")
#convert rasterStack to dataframe
stackPoints <- as.data.frame(rasterToPoints(stack))
#format x and x coordinates
colX <- formatCoordinate(stackPoints$x)
colY <- formatCoordinate(stackPoints$y)
#combine formatted x and y coordinates to compose a unique key
pK <- paste0(colX, colY )
stackPoints["key"] <- pK
col_idx <- grep("key", names(stackPoints))
stackPoints <- stackPoints[, c(col_idx, (1:ncol(stackPoints))[-col_idx])]
#write results to a csv file
write.table(stackPoints, "r.csv", row.names=F, sep=";", dec=",", append=F)

上面的代码适用于小型栅格，但对于大型栅格，我无法将堆栈加载到 RAM。 有没有办法将我的代码转换为使用并行处理？即使用多核读取光栅并写入 csv，而无需将光栅加载到 RAM（Mac OSX 10.11 和 Ubuntu 14.04，每个 8 核）。 最好的，

Answer 1

首先，您想弄清楚如何在单个线程上编写循环，因为从for() 移动到foreach() 将非常简单。我不熟悉 RasterStack 对象，但它看起来像它的层可以用nlayers(x) 计数并且可以用x[[i]] 提取。

所以首先我会编写和调试类似的东西：

for(i in 1:nlayers(stack)){
  #convert layer of rasterStack to dataframe
  layer_pts <- as.data.frame(rasterToPoints(stack[[i]]))

  #write layer_pts to a csv file
}

那么foreach() 就很简单了。请记住，您需要使用 raster 包启动每个线程。对于更快的合并，我推荐data.table。

library(foreach)
library(doMC)
library(data.table)
registerDoMC(detectCores() - 2) # for me this is 40 - 2 = 38
layer_list <- 
  foreach(i = 1:nlayers(stack), .packages = c('raster', 'data.table') ) %dopar% {
    #convert layer of rasterStack to data.table
    layer_pts <- as.data.table(rasterToPoints(stack[[i]]))
    setkey(layer_pts, x, y) # data.table can key on x and y, no synthetic key needed
    layer_pts
  }

tbl_out <- Reduce(merge, layer_list) # uses keys from setkey

# if you wanted the "key" column (but not essential)
tbl_out[, key:= paste0( formatCoordinate(x), formatCoordinate(y) ) ]

write.csv(tbl_out, 'r.csv')

请注意，如果内存不足，您可能必须减少使用的内核数量。例如，registerDoMC(4) 基于反复试验。