根据 RAM 和速度在 r 中的平衡优化栅格到数据帧函数

Question

问题

我正在尝试使用一个模型，它需要我将非常大的 rastersStacks（大约 1000 万个单元格）转换为 data.frame，我在共享服务器上执行此操作，因此，我我正在尝试优化以减少使用的 RAM，并希望不会大幅降低速度。为此，我编写了 2 个函数，但都没有成功。包含我尝试结果的 RPUBS 在此 link, and a github with the rmd for that is here 中，包括用于 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 个参数：

Raster:你要分块的堆栈
ppside:你最终得到的水平和垂直瓦片的数量，最终的瓦片数量为ppside*ppside，所以如果ppside为3，你将有9瓷砖
nclus: 用于并行化和加速处理的集群数。

在此函数结束时，您将拥有 ppside*ppside 个图块，保存在您的工作目录中，所有名称都称为 SplitRasN.tif，其中 N 是图块的编号。举个例子，我们将使用栅格包中可用的生物气候变量：

Bios <- getData('worldclim', var='bio', res=10)

我将对此进行测试，将其转换为不同数量的图块，如下图所示

从栅格到切片然后从切片到数据框的转换（示例）

所以首先我们将使用 SplitRas 函数使用 4 个核心获得 16 个图块：

SplitRas(Raster = Bios, ppside = 4, nclus = 4)

这将 return 以下文件：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")

Answer 1

您可以使用 rasterToPoints。请注意，输出中省略了全部为 NA 的任何行。还值得指出的是，您可以通过 rasterOptions(maxmemory = XXX).

控制使用多少 RAM

x = as.data.frame(rasterToPoints(Bios))
head(x)
#           x        y bio1 bio2 bio3  bio4 bio5 bio6 bio7 bio8 bio9 bio10 bio11 bio12 bio13 bio14 bio15 bio16 bio17 bio18 bio19
# 1 -37.91667 83.58333 -174   67   17 11862   37 -356  393  -31 -307    -7  -307   144    22     7    38    59    24    50    24
# 2 -37.75000 83.58333 -174   67   17 11870   37 -355  392  -30 -219    -7  -307   143    22     7    42    59    23    50    24
# 3 -36.91667 83.58333 -172   68   17 11872   39 -354  393  -29 -217    -5  -305   136    22     6    42    57    22    49    23
# 4 -36.75000 83.58333 -173   68   17 11887   39 -354  393  -29 -217    -5  -306   136    22     6    42    57    22    49    23
# 5 -36.58333 83.58333 -173   68   17 11877   39 -354  393  -29 -217    -6  -306   136    22     6    42    57    22    49    23
# 6 -36.41667 83.58333 -173   68   17 11879   38 -354  392  -29 -217    -6  -306   137    22     6    41    57    22    49    24

Answer 2

I am trying to use a model which requires me to transform, very large rastersStacks (around 10 million cells) to a data.frame,

相反，我会使用 raster::predict 或 terra::predict；也许运行这些具有并行化。

terra 有一个 makeTiles 方法可能有用。

根据 RAM 和速度在 r 中的平衡优化栅格到数据帧函数

Optimization of raster to data frame function in terms of balance of RAM and speed in r

ram

r

raster

r-raster

profvis

问题

第一个函数

平铺

从栅格到切片然后从切片到数据框的转换（示例）

内存基准测试（我试过的）

一些额外的东西