问题
我正在尝试使用一个模型，该模型需要我将非常大的 rasterStacks（大约 1000 万个单元）转换为 data.frame，我在共享服务器上执行此操作，因此，我正在尝试优化以减少使用了 RAM，希望不会大大降低速度。为此，我编写了 2 个函数，但没有成功。带有我尝试结果的 RPUBS 位于此链接中，带有 rmd 的 github 位于此处，包括用于 profvis 分析的 rds 文件。
第一个功能
首先让我们加载我们需要的包：
# For spaital manipulation
library(raster)
# For benchmarking speed and memory
library(profvis)
# To parallelize operations
library(doParallel)
# To parallelize operations
library(foreach)
# For combination and looping
library(purrr)
# Data wranggling
library(dplyr)
library(data.table)
平铺
减少 RAM 使用的主要方法不是处理一个大栅格，而是将其转换为平铺栅格，为此我开发了以下功能：
SplitRas <- function(Raster,ppside, nclus = 1){
TempRast <- paste0(getwd(), "/Temp")
h <- ceiling(ncol(Raster)/ppside)
v <- ceiling(nrow(Raster)/ppside)
agg <- aggregate(Raster,fact=c(h,v))
agg[] <- 1:ncell(agg)
agg_poly <- rasterToPolygons(agg)
names(agg_poly) <- "polis"
r_list <- list()
if(nclus == 1){
for(i in 1:ncell(agg)){
dir.create(TempRast)
rasterOptions(tmpdir=TempRast)
e1 <- extent(agg_poly[agg_poly$polis==i,])
r_list[[i]] <- crop(Raster,e1)
if((freq(r_list[[i]], value=NA)/ncell(r_list[[i]])) != 1){
writeRaster(r_list[[i]],filename=paste("SplitRas",i,sep=""),
format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
unlink(TempRast, recursive = T, force = T)
}
} else if(nclus != 1){
cl <- parallel::makeCluster(nclus)
doParallel::registerDoParallel(cl)
r_list <- foreach(i = 1:ncell(agg), .packages = c("raster")) %dopar% {
dir.create(TempRast)
rasterOptions(tmpdir=TempRast)
e1 <- extent(agg_poly[agg_poly$polis==i,])
Temp <- crop(Raster,e1)
if((raster::freq(Temp, value=NA)/ncell(Temp)) != 1){
writeRaster(Temp,filename=paste("SplitRas",i,sep=""),
format="GTiff",datatype="FLT4S",overwrite=TRUE)
}
unlink(TempRast, recursive = T, force = T)
Temp
}
parallel::stopCluster(cl)
}
}
我的想法是，如果您分别处理每个图块，则可以一个一个地转换为数据帧，并取出带有 NA 的行，从而减少 RAM 使用量。
该函数有 3 个参数：
- 光栅：你想分成瓷砖的堆栈
- ppside：你最终会得到的水平和垂直瓷砖的数量，最终的瓷砖数量将是 ppside*ppside，所以如果 ppside 为 3，你将有 9 个瓷砖
- nclus：您将用于并行化和加速处理的集群数量。
在此函数结束时，您将拥有ppside*ppside多个切片，保存在您的工作目录中，所有这些都称为SplitRasN.tif其中 N 是切片的编号。作为示例，我们将使用 raster 包中可用的生物气候变量：
Bios <- getData('worldclim', var='bio', res=10)
我将对此进行测试，将其转换为不同数量的图块，如下图所示
从栅格转换为切片，然后从切片转换为数据框（示例）
所以首先我们将使用SplitRas函数来获得使用 4 个核心的 16 个图块：
SplitRas(Raster = Bios, ppside = 4, nclus = 4)
这将返回以下文件：r list.files(pattern = "SplitRas")
为了将这些图块放入一个包含所有非 NA 单元格的数据帧中，我们需要一个图块列表，我们可以通过以下方式获得：
Files <- list.files(pattern = "SplitRas", full.names = T)
我们可以在以下函数中使用它：
SplitsToDataFrame <- function(Splits, ncores = 1){
TempRast <- paste0(getwd(), "/Temp")
if(ncores == 1){
Temps <- list()
for(i in 1:length(Splits)){
dir.create(TempRast)
rasterOptions(tmpdir=TempRast)
Temp <- stack(Splits[i])
Temp <- as.data.frame(Temp, row.names = NULL, col.names = NULL, xy =TRUE)
colnames(Temp)[3:ncol(Temp)] <- paste0("Var", 1:length(3:ncol(Temp)))
Temps[[i]] <- Temp[complete.cases(Temp),]
gc()
unlink(TempRast, recursive = T, force = T)
message(i)
}
Temps <- data.table::rbindlist(Temps)
} else if(ncores > 1){
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
Temps <- foreach(i = 1:length(Splits), .packages = c("raster", "data.table")) %dopar%{
dir.create(TempRast)
rasterOptions(tmpdir=TempRast)
Temp <- stack(Splits[i])
Temp <- as.data.frame(Temp, row.names = NULL, col.names = NULL, xy =TRUE)
colnames(Temp)[3:ncol(Temp)] <- paste0("Var", 1:length(3:ncol(Temp)))
gc()
unlink(TempRast, recursive = T, force = T)
Temp[complete.cases(Temp),]
}
Temps <- data.table::rbindlist(Temps)
parallel::stopCluster(cl)
}
return(Temps)
}
其中Splits是具有图块路径的字符向量，ncores是用于并行化的核心数。这将导致具有非空单元格的数据框。
DF <- SplitsToDataFrame(Splits = Files, ncores = 4)
内存基准测试（我尝试过的）
首先，我为 1、2、4、8、10 和 12 ppsides 生成了 Tiles
sides <- c(1,2,4,8,10, 12)
Home <- getwd()
AllFiles <- list()
for(i in 1:length(sides)){
dir.create(path = paste0(Home, "/Sides_", sides[i]))
setwd(paste0(Home, "/Sides_", sides[i]))
SplitRas(Raster = Bios, ppside = sides[i], nclus = ifelse(sides[i] < 4, sides[i], 4))
AllFiles[[i]] <- list.files(pattern = "SplitRas", full.names = T) %>% stringr::str_replace_all("\\./", paste0(getwd(), "/"))
}
setwd(Home)
然后仅使用函数的顺序选项来分析函数
library(profvis)
P <- profvis({
P1 <- SplitsToDataFrame(Splits = AllFiles[[1]])
P2 <- SplitsToDataFrame(Splits = AllFiles[[2]])
P3 <- SplitsToDataFrame(Splits = AllFiles[[3]])
P4 <- SplitsToDataFrame(Splits = AllFiles[[4]])
P5 <- SplitsToDataFrame(Splits = AllFiles[[5]])
})
P
htmlwidgets::saveWidget(P, "profile.html")
saveRDS(P, "P.rds")
然而，如下图所示（更多细节可以在上面链接的 Rpubs 中找到），RAM 或多或少与以前相同，但时间变长了（最后一部分是预期的）。
.
一些额外的东西
最后，当我尝试并行对 RAM 使用进行基准测试时，似乎 profvis 没有捕捉到这一点。关于如何检查的任何想法？
library(profvis)
PPar <- profvis({
P1 <- SplitsToDataFrame(Splits = AllFiles[[3]], ncores = 1)
P2 <- SplitsToDataFrame(Splits = AllFiles[[3]], ncores = 2)
P3 <- SplitsToDataFrame(Splits = AllFiles[[3]], ncores = 4)
P4 <- SplitsToDataFrame(Splits = AllFiles[[3]], ncores = 7)
})
PPar
htmlwidgets::saveWidget(PPar, "profileParallel.html")
saveRDS(PPar, "PPar.rds")
