GRASS GIS 中的 r.mask 和 readRAST6

生成上篇记录里那种印花桌布的原因已经查明，修改后的脚本正在跑，应该最多只有些小问题了。

我早该想到的，就是我以为在 GRASS 中对图像进行 mask 操作以后，用 as.vector(data.frame(readRAST6())) 这样的操作得到的图像仍然是全幅的，缺少数值的地方用 NA 填充。事实上这个操作去掉了所有标记为 NA 的像素，留给我的只有少数几个像素而已。R 在这种情况下默认用重复较短的数组/矩阵来填充较长的，于是我就得到了漂亮的印花桌布-_-

事实上 readRAST6() 本身返回的结果就是类似于 data frame 的结构，其中第一个变量是我们关心的属性，假设原来的栅格图像在 GRASS 中的名字叫 foo，导入 R 后的结构名称为 bar，则 bar 中这个变量的名字可以用 names(foo) 获得，可能是类似 "Xfoo" 这样的形式。在 R 下用 bar[["Xfoo"]] 这样的形式引用就可以了（bar$Xfoo 大概也行）。

早上就写完了脚本，因为这个问题 debug 到半夜，当然主要还是对 R 不熟悉，其间也发现了一些脚本其他的问题。最近看一本国人写的 R 教材，有些我使用中知道的东西书上也没有提到。莫非 R 已经庞大到无法系统全面地描述其所有的功能了么？spgrass6 这个函数包的文档相当简陋，更让我有 R 的发展超过了文档编写的速度的感觉，这未必是件好事……

这个脚本总共要计算 540 * 7 个图像，等它算完估计天亮了……我先去见周公，明天看结果，如果图好看的话补贴图。

将 R 中的向量作为空间数据导入 GRASS GIS

GRASS GIS 6.x 中的栅格和矢量数据可以通过 spgrass6 这个 R 的扩展包导入 R 进行统计和其他计算。用 readRAST6() 导入的栅格数据是一个称为 SpatialGridDataFrame 的结构，含有栅格位置、尺寸和投影等信息，栅格的属性数据则按照先从左到右，再从上到下的扫描顺序排成一个列向量，和对应的坐标位置一起储存在结构中。为了方便计算，这个向量可以用类似 as.vector(data.frame(readRAST6())[,1]) 这样的方式提取为普通的数值向量（is.numeric()=TRUE）。但运算完毕后要再将其导出到 GRASS GIS 中，由于缺少空间信息，略微有些麻烦。

这里提供了一些有用的讨论来解决这个问题。虽然向量不能直接导入 GRASS，但是向量可以赋值给一个 SpatialGridDataFrame，例如这样：

ensoeffarea<-array(0,dim<-c(13640,3))
ensoeffarea[,1]<-ndvimaxvalues-ndvielvalues
ensoeffarea[,2]<-data.frame(readRAST6(as.character(filename[1,1])))[,2]
ensoeffarea[,3]<-data.frame(readRAST6(as.character(filename[1,1])))[,3]
temp<-readRAST6(as.character(filename[1,1]))
temp$temp<-ensoeffarea[,1]
writeRAST6(temp,"enso-strike-area","temp")

要注意用于赋值的向量必须是列向量。而且临时使用的这个 SpatialGridDataFrame 的各种空间信息必须与向量所要求的一致，特别是各像素或特征对象数值与像素位置/对象 id 的对应关系不应有变化。否则就只能采用比较费劲的先导出成文本文件再修改导入 GRASS 的方法了。