rgdal高效读取大型多波段光栅

我正在使用rgdal包在R中编写一个图像分类脚本。所讨论的光栅是一个PCIDSK文件,有28个通道:一个训练数据通道、一个验证数据通道和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所建议的,使用“光栅”软件包要快得多;但是据我所知,就收集训练数据而言,它并不能解决内存问题,因为它最终需要在一个数据帧中,在内存中:

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)

因此,虽然这要快得多(需要的代码也更少),但它仍然不能解决没有足够的空闲内存来创建分类器的问题。。。是否有一些“光栅”包功能,我没有发现,可以做到这一点?谢谢。。。