小功告成

在熬了两夜，整整两天连登陆学校邮箱都顾不上的不断 scripting 和 debug 之后，我总算完成了研究地区的一部分净初级生产力分布图。

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 个图像，等它算完估计天亮了……我先去见周公，明天看结果，如果图好看的话补贴图。

Pattern by Error

用 bash 和 R 联合编了一个比较复杂的脚本，计算研究区域的 FAPAR （Fraction of Absorbed Photosynthetically Active Radiation）。但是在某一步出错了，产生了这样的东西……

将 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 的方法了。

A flavour of AVHRR

In the past few days I have been trying to find suitable software and process for AVHRR data processing. The METOP AVHRR format looks like an enigma. It should be compatible to AVHRR, but despite that I tried hard to find a tool to unpack the file, and even tried to write a Fortran program myself, I couldn't find anything that decodes the file downloaded from [CLASS](http://www.class.noaa.gov/).

What made the search not totally fruitless is a image processing tool, [BEAM](http://www.brockmann-consult.de/beam/), which works nicely with AVHRR and a number of other data files. I can use it to open AVHRR 16-bit and 8-bit Level 1b file (10-bit not tested yet), do some simple preprocessing and statistics, then export the file as GeoTIFF or other formats, which can be read by GRASS GIS + GDAL. It is worth noting that BEAM is free as in free speech software published under GPL **v3** and can run on any platform with Java Runtime Environment.

With the help of BEAM and GRASS I managed to produce my first AVHRR-derived NDVI image this afternoon. Seems I am finally on track. In the coming days there are a number of things that I can look at:

- AVHRR Pathfinder data, their advantages and limits, suggested by one of my supervisor;
- Snow and ice cover mapping with AVHRR;
- Time series data and change detection;
- MODIS data and its interoperability with AVHRR.

自由的遥感图像处理解决方案？

因为工作的需要，最近一段时间的工作重点是搜集和处理遥感图像。具体要怎么处理还没有决定，只知道大概要做的东西包括几何校正，大气校正，geoencoding/registration，后期肯定会做 change detection。这些在 ERDAS Imagine 里面是相对简单的，基本都可以点点鼠标解决。但是目前学校的 IT staff 还没有给我的电脑装这套软件，为了立即开展工作我准备找一找自由/开源的替代方案。

在 [http://opensourcegis.org/](http://opensourcegis.org/) 上给出了一个比较全面的与 GIS/RS 有关的自由/开源项目的列表，虽然其中成气候的项目还很少，而且像 Imagemagick 和 Free Type 这样的通用项目都列上去了，不过这方面的开发工作确实比我想象的更多。在 [GDAL](http://www.gdal.org/) 出来后出现了很多 GIS 项目，但地学图像处理方面的软件相对还是很缺乏。

一番比较后我决定使用 [GRASS GIS](http://grass.itc.it/intro/general.php) 和 [nip2](http://www.vips.ecs.soton.ac.uk/index.php?title=Nip2) 来完成工作。GRASS GIS 是目前自由/开源 GIS 中最成熟的项目（应该不用加“之一”了吧）。以前我用它做过一个简单的模型输出可视化，但没注意到它还有处理遥感图像的功能。看了它的手册和一些搜索之后，结论是它可以用来做几何与大气校正、合成和分类，虽然分类的算法可能没有商业软件包那么多样，不过理论上说可以自己写个实现，已经有人做出了[神经网络分类方法的实现](http://www.ncgia.ucsb.edu/conf/SANTA_FE_CD-ROM/sf_papers/muttiah_ranjan/muttiah.html)。

nip2 是 [VIPS](http://www.vips.ecs.soton.ac.uk/index.php?title=VIPS) 库的一个前端，号称是介于 Photoshop 与 Excel 之间的图像处理软件。其特点是可以用控件和脚本结合对图像方便地进行复杂的处理。虽然它的设计似乎并没有特别考虑地学应用的需要，但它有一个非常吸引人的特性，就是对大图像的处理速度非常快，而且可以处理大于内存容量的图像。我准备用它来处理观测角度与读数关系不大的图像，也可能用它来做分割和组合。

如果两周后我觉得这两个软件用得比较顺手，可能会一直用它们完成论文。商业软件包确实易用，能够节省一些学习和适应的时间，但我觉得用自由/开源的软件能让自己对 under the hood 的东西了解更多一些，而且现在掌握了 GRASS GIS 以后还可以很快上手 ERDAS Imagine，反之则不然。