rgdal 有效地读取大型多波段栅格答案

【问题标题】：rgdal efficiently reading large multiband rastersrgdal 有效地读取大型多波段栅格
【发布时间】：2010-11-15 16:28:17
【问题描述】：

我正在使用 rgdal 包在 R 中编写图像分类脚本。有问题的栅格是一个具有 28 个通道的 PCIDSK 文件：一个训练数据通道、一个验证数据通道和 26 个光谱数据通道。目标是填充一个数据帧，其中包含训练数据通道中不为零的每个像素的值，以及 26 个波段中的相关光谱值。

在 Python/Numpy 中，我可以轻松地将整个图像的所有波段导入到多维数组中，但是，由于内存限制，R 中唯一的选择似乎是逐块导入此数据，这非常慢：

library(rgdal)

raster = "image.pix"
training_band = 2
validation_band = 1
BlockWidth = 500
BlockHeight = 500

# Get some metadata about the whole raster
myinfo = GDALinfo(raster)
ysize = myinfo[[1]]
xsize = myinfo[[2]]
numbands = myinfo[[3]]


# Iterate through the image in blocks and retrieve the training data
column = 0
training_data = NULL
while(column < xsize){
 if(column + BlockWidth > xsize){
  BlockWidth = xsize - column
 }
 row = 0
 while(row < ysize){
  if(row + BlockHeight > ysize){
   BlockHeight = ysize - row
  }

  # Do stuff here
  myblock = readGDAL(raster, region.dim = c(BlockHeight,BlockWidth), offset = c(row, column), band = c(training_band,3:numbands), silent = TRUE)
  blockdata = matrix(NA, dim(myblock)[1], dim(myblock)[2])
  for(i in 1:(dim(myblock)[2])){
   bandname = paste("myblock", names(myblock)[i], sep="$")
   blockdata[,i]= as.matrix(eval(parse(text=bandname)))
  }
  blockdata = as.data.frame(blockdata)
  blockdata = subset(blockdata, blockdata[,1] > 0)
  if (dim(blockdata)[1] > 0){
   training_data = rbind(training_data, blockdata)
  }

  row = row + BlockHeight
 }
 column = column + BlockWidth
}
remove(blockdata, myblock, BlockHeight, BlockWidth, row, column)

有没有更快/更好的方法来做同样的事情而不会耗尽内存？

收集此训练数据后的下一步是创建分类器（randomForest 包），该分类器也需要大量内存，具体取决于请求的树数。这就引出了我的第二个问题，即考虑到训练数据已经占用的内存量，不可能创建一个由 500 棵树组成的森林：

myformula = formula(paste("as.factor(V1) ~ V3:V", dim(training_data)[2], sep=""))
r_tree = randomForest(formula = myformula, data = training_data, ntree = 500, keep.forest=TRUE)

有没有办法分配更多的内存？我错过了什么吗？谢谢...

[编辑] 正如 Jan 所建议的，使用“raster”包要快得多；但是据我所知，就收集训练数据而言，它并不能解决内存问题，因为它最终需要在数据帧中，在内存中：

library(raster)
library(randomForest)

# Set some user variables
fn = "image.pix"
training_band = 2
validation_band = 1

# Get the training data
myraster = stack(fn)
training_class = subset(myraster, training_band)
training_class[training_class == 0] = NA
training_class = Which(training_class != 0, cells=TRUE)
training_data = extract(myraster, training_class)
training_data = as.data.frame(training_data)

因此，虽然这要快得多（并且需要更少的代码），但它仍然不能解决没有足够的可用内存来创建分类器的问题......是否有一些我没有发现的“光栅”包函数能做到这一点吗？谢谢...

【问题讨论】：

标签： r raster gdal

【解决方案1】：

查看光栅包。 Raster 包为 Rgdal 提供了一个方便的包装器，无需将其加载到内存中。

http://raster.r-forge.r-project.org/

希望对您有所帮助。

“光栅”包处理基本的空间栅格（网格）数据访问和操纵。它定义了栅格课程；可以处理非常大文件（存储在磁盘上）；并包括标准栅格函数，例如覆盖、聚合和合并。

“光栅”包的目的是提供易于使用的功能光栅操作和分析。这些包括高级功能如叠加、合并、聚合、投影，重采样，距离，多边形到栅格的转换。全部这些功能适用于非常大的无法加载的栅格数据集进入记忆。此外，包提供较低级别的功能，例如逐行读写（到通过 rgdal 的多种格式）用于构建其他功能。

通过使用 Raster 包，您可以避免在使用 randomForest 之前填满内存。

[编辑] 如果您可以在子样本（大小

【讨论】：

感谢您的回复。我看了一下“光栅”包，确实它更合适。但是，在收集训练数据方面，我仍然遇到同样的问题：需要将其加载到数据框中才能创建分类器。请参阅我的原始问题以获取我的新代码。谢谢，本杰明。
我会按照你的建议单独生成树。

【解决方案2】：

我认为这里的关键是：“一个数据帧，其中包含训练数据通道中不为零的每个像素的值”。如果生成的 data.frame 足够小以保存在内存中，您可以通过仅读取该波段来确定这一点，然后只修剪那些非零值，然后尝试创建一个具有那么多行和总数的 data.frame你想要的列。

你能运行这个吗？

 training_band = 2 

 df = readGDAL("image.pix",  band = training_band) 
 df = as.data.frame(df[!df[,1] == 0, ])

然后，您可以通过分别读取每个波段并修剪训练波段来一一填充 data.frame 的列。

如果那个 data.frame 太大，那么你就卡住了——我不知道 randomForest 是否可以使用“ff”中的内存映射数据对象，但它可能值得一试。

编辑：一些示例代码，并注意 raster 为您提供内存映射访问，但问题是 randomForest 是否可以使用内存映射数据结构。您可能只能读取您需要的数据，一次一个波段 - 您可能希望首先尝试构建完整的 data.frame，而不是附加列。

此外，如果您可以从一开始就生成完整的 data.frame，那么您就会知道它是否应该工作。通过 rbind() 处理代码，您需要越来越大的连续内存块，这可能是可以避免的。

【讨论】：