首页 文章

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

提问于
浏览
2

我正在使用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所建议的,使用“光栅”包的速度要快得多;但据我所知,就收集训练数据而言,它并没有解决内存问题,因为它最终需要在数据帧中,在内存中:

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)

因此虽然速度更快(并且代码更少),但它仍然无法解决没有足够的可用内存来创建分类器的问题...是否有一些我没有找到的“光栅”包函数可以实现这一点?谢谢...

2 回答

  • 3

    查看Raster包 . Raster包为Rgdal提供了一个方便的包装器,而无需将其加载到内存中 .

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

    希望这有帮助 .

    'raster'包处理基本的空间栅格(网格)数据访问和操作 . 它定义了栅格类;可以处理非常大的文件(存储在磁盘上);并包括标准栅格功能,如叠加,聚合和合并 . “光栅”包的目的是为光栅操作和分析提供易于使用的功能 . 这些包括高级功能,如叠加,合并,聚合,投影,重采样,距离,多边形到光栅转换 . 所有这些函数都适用于无法加载到内存中的超大型栅格数据集 . 此外,该软件包还提供较低级别的功能,例如逐行读取和写入(通过rgdal的多种格式)以构建其他功能 .

    通过使用Raster包,您可以避免在使用randomForest之前填充内存 .

    [编辑]为了解决randomForest的内存问题,如果您可以在子样本(大小为<< n)而不是引导样本(大小为n)上学习随机林中的单个树,则可能会有所帮助 .

  • 0

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

    你能跑吗?

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

    然后,您可以通过分别读取每个波段并对训练波段进行修整来逐个填充data.frame的列 .

    如果那个data.frame太大,那你就被卡住了 - 我不知道randomForest是否可以在“ff”中使用内存映射数据对象,但它可能值得尝试 .

    编辑:一些示例代码,并注意栅格为您提供内存映射访问,但问题是randomForest是否可以使用内存映射数据结构 . 您可以只读取所需的数据,一次只读取一个 - 您可能希望首先构建完整的data.frame,而不是追加列 .

    此外,如果您可以从一开始就生成完整的data.frame,那么您将知道它是否应该正常工作 . 通过rbind()在您的代码执行过程中,您需要越来越大的连续内存块,这可能是可以避免的 .

相关问题