SRTP结题报告

西南交通大学第五期大学生科研训练计划（SRTP）

遥感与气象观测数据耦合的地表干

旱状态监测

结题报告

2010年4月至2011年4月

目录

目录_______________________________________________________________ 1

1绪论 ____________________________________________________________ 2

1. 1项目背景 _________________________________________________________________ 2

1. 2项目介绍 _________________________________________________________________ 3

1. 3 开源库GDAL与HDF文件格式介绍________________________________________ 4

2数据源及GDAL库的配置 __________________________________________ 7

2. 1 数据源__________________________________________________________________ 7

2. 2 GDAL库的配置 ___________________________________________________________ 9

3数据处理原理及实现 _____________________________________________ 10

3. 1 遥感数据处理___________________________________________________________ 10

3.2 气象数据处理与降水空间插值法比较 _____________________________________ 20

3. 3 指数耦合 _______________________________________________________________ 26

3. 4 程序设计 _______________________________________________________________ 27 4实验结果分析 ___________________________________________________ 27 5不足与期望 _____________________________________________________ 30 6项目感想 _______________________________________________________ 31参考文献__________________________________________________________ 32附录1 程序主要源代码 ______________________________________________ 34

1绪论

1.1项目背景

我国是一个旱灾非常严重的国家，旱灾给农业、农村和农民造成了巨大的损失。据统计，1950—1990年间，我国共有11年发生了重特大干旱，发生频次为26.83％，因干旱造成粮食损失占粮食总产量的4.02％。而1991—2009年间，我国共有8年发生重特大干旱，因干旱造成粮食损失占粮食总产量的6.09％。近年来，我国年年有干旱，平均不到3年发生一次重特大旱灾，尤其经常发生区域性特大旱灾。同时我国旱灾分布面积广，过去，我国旱灾高发的区域主要在干旱缺水的北方地区，特别是西北地区。近几年，在传统的北方旱区旱情加重的同时，南方和东部多雨区旱情也在扩展和加重，目前旱灾范围已遍及全国。与此同时，旱灾影响范围已由传统的农业扩展到工业、城市、生态等领域，工农业争水、城乡争水、超采地下水和挤占生态用水现象越来越严重。因此实现对旱灾的监测，有利于实施抗旱措施，同时合理分配水资源，节约用水，对促进农业生产、保障粮食安全和区域可持续发展具有重要的现实意义1。

传统的干旱监测采用气象数据，气象数据由气象站观测得到，通过建立区域气象站网同时辅以水文、社会和经济等数据实施干旱监测。这种方法应用十分广泛，在单点精度很高，但其性质是以点带面来监测，这样很难给出不同干旱区域的分界线，同时整体精度也相对较低，难以适应大范围的干旱监测。

随着卫星遥感技术的迅速发展，国内外已在借助遥感手段进行大范围的干旱监测方面开展了大量的研究和实际应用。基于遥感手段的干旱监测能频繁和持久的提供地表特征的面状信息，具有宏观，高时效，经济等特点，且点点俱到，适应大范围的干旱监测。当然，由于气象卫星运行特点和遥感传感器本身性能的限制，该方法也存在一定的缺点和问题，比如目前相应的基于遥感数据的一些反演模型不够成熟，并且其受大气影响较大，在很多有云的情况下，数据不能使用，其监测的结果在小范围内也不如传统的气象监测精确。

1胡其峰.《科学理性对待西南旱灾减灾需立足长远》. 光明日报.2010.4

针对单一遥感数据在进行地表干旱状态监测时的不足，考虑降水对土壤水分的阶段性影响，我们应该把两种方法结合在一起考虑，形成地表温度、植被状态、空间降水等参数耦合的干旱监测模型。这样能提高干旱监测的准确性和及时性。

1.2项目介绍

项目基于对地观测MODIS数据与气象站点实测资料，建立起多源数据耦合下的地表干旱监测方法与相关指标体系，在开源GDAL库的基础上，开发可用于大范围区域地表干旱状态监测的应用系统。

项目中涉及四个方面的内容：MODIS预处理算法程序设计、气象观测资料处理、热红外地表温度反演模型程序设计、多元数据耦合模型参数化方案及设计。其中，MODIS预处理算法程序设计部分包括HDF格式数据读取、MODIS几何校正与辐射标定处理。气象观测资料处理包括气象观测站点降水资料读取、空间插值与栅格化。热红外地表温度反演程序设计主要实现利用分裂窗法的地表温度反演模型。耦合模型部分则包括了地表温度参数与降水参数的耦合模型参数计算、模型程序实现。

本项目主要分为三个部分：遥感数据处理，气象数据处理以及干旱监测指数的耦合。

遥感数据处理：利用开源库GDAL提取MODIS 1B的数据，在此基础上利用MODIS数据中的相应波段反演地表温度（LST）和归一化植被指数（NDVI），得到遥感干旱监测因子作物供水指数（VSWI）。

气象数据处理：利用ArcGIS的二次开发平台ArcGIS Engine实现气象数据降雨量的多种空间插值算法。在比较不同空间插值算法误差的基础上，综合多年数据选取最佳的空间插值算法进行降雨量的插值运算，并计算出气象干旱监测因子综合降水距平指数（MSRI）。

干旱监测指数的耦合：建立作物供水指数和综合降水距平指数的耦合模型，得到农业旱情监测指数，并划分干旱等级。

项目流程图如图1所示：

图1-1项目流程图

1. 3 开源库GDAL 与HDF 文件格式介绍

1. 3. 1 GDAL 库介绍

在MODIS 数据的处理中，我们没有选用HDF Group 开发的HDF 函数库，而是选用了十分流行的开源库GDAL 。因为开源库GDAL 在处理栅格图像格式有很大优势，其不仅读取数据较快，而且能够把结果转换成任一其支持的栅格格式。

1998年末，加拿大的Frank Warmerdam 开始了GDAL(Geospatial Data

Abstraction Library)项目的编写工作。该项目得到了许多个人和团队的支持。GDAL 是一个在X/MIT 许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。OGR 是GDAL 项目的一个分支，功能和GDAL 类似，只不过他提供对矢量数据的支持。因此GDAL 能同时提供栅格和矢量数据的操作。有很多著名的GIS 类产品包括ESRI 的ARCGIS ，Google Earth 和跨平台的GRASS GIS 系统都使用了GDAL/OGR 开源库的若干模块。利用GDAL/OGR 库，可以使基于Linux 的地理空间数据管理系统提供对矢量和栅格文件数据的支持。

遥感数据->MODIS1B 数据 地面观测数据->降雨量

云检测

NDVI LST

VSWI ＝NDVI/LST

单日数据内插

旬合成

当旬降水量距平

综合降水距平指数

旬合成

综合旱情监测指数

几何校正

GDAL提供对多种栅格数据的支持，包括Arc/Info ASCII Grid(asc)，GeoTiff (tiff)，Erdas Imagine Images(img)，ASCII DEM(dem) 等格式。提供读取、写入、转换、处理各种栅格数据格式（有些特定的格式对一些操作如写入等不支持）操作。GDAL使用抽象数据模型(Abstract Datamodel)来解析它所支持的数据格式，抽象数据模型包括数据集(Dataset)，坐标系统，仿射地理坐标转换(Affine GeoTransform)，大地控制点(GCPs)，元数据(Metadata)，栅格波段(Raster Band)，颜色表(ColorTable)，子数据集域(Subdatasets Domain)，图像结构域

(Image_StructureDomain)，XML域(XML:Domains)。

GDAL库采用ANSI C和C++语言编写，能通过所有c/c++编译器的编译。其核心类框架如图2所示

GDALMajorObject

GDALDataset GDALDrive GDALDriverManager GDALRasterBand

图1-2 GDAL核心类框架

其中：

GDALMajorObject类：带有元数据的对象。

GDALDataset类：通常是从一个文件中提取的相关栅格波段数据的集合，以及其相关元数据。

GDALDriver类：文件格式驱动类。

GDALDriverManager类：文件格式驱动管理类。

GDALRasterBand类：栅格波段数据类。

1. 3. 2 HDF文件格式

HDF（Hierarchical data format）格式是美国国家高级计算机应用中心为了满足各种领域研究需求而研制的一种高效存储和分发科学数据的层次数据格式，主

要用来存储由不同计算机平台来产生的各种类型科学数据，便于在不同的计算机平台上扩展。HDF 被设计为：

自述性：对于一个HDF 文件里的每一个数据对象，有关于该数据的综合信息。 通用性：许多数据类型都可以被嵌入在一个HDF 文件里。

灵活性：HDF 允许用户把相关的数据对象组合在一起，放到一个分层结构中，向数据对象添加描述和标签。它还允许用户把科学数据放到多个HDF 文件里。

扩展性：HDF 极易容纳将来新增加的数据模式，容易与其他标准格式兼容。 跨平台性：HDF 是一个与平台无关的文件格式。HDF 文件无需任何转换就可以在不同平台上使用。

目前HDF 有两个版本，即HDF4和HDF5，两个版本的数据结构模型变化不大，但其具体实现方法变化很大，因此这两个版本是不兼容的。由于GDAL 库采用一个抽象的数据模型来支持遥感栅格图像，因此使用GDAL 读取这两个版本的HDF 数据的方法基本相同。一个HDF 文件包括一个文件头，至少一个数据描述块以及若干个数据元素（图1-3）。HDF 文件像一本带目录的书，数据描述块可以看成书的目录，而书的每一章内容则是HDF 文件的数据元素。

图1-3 HDF 文件格式

……

数据元素

HDF 文件头

……

空数据描述符

数据描述符

空数据描述符

空数据描述符

HDF 文件 数据描述块

本项目采用的MODIS 1B数据采用HDF文件格式储存。图1-4 是通过HDFView软件浏览器查看某MODIS 1km分辨率的文件列表示意图。

波段：15*2040*1354

波段：2*2040*1354

波段：5*2040*1354

卫星成像参数：

408*271

图1-4 HDF（MODIS）文件列表示意图

2数据源及GDAL库的配置

2.1数据源

2. 1. 1遥感数据源

中等分辨率成像光谱仪MODIS数据采用HDF(Hierarchical data format)数据格式存储元数据，提供36个波段的地球综合信息，分布在0.4~14μm的电磁波谱范围内，对开展自然灾害与生态环境监测、全球环境和气候变化研究以及进行全球变化的综合性研究等有重要意义。MODIS仪器的地面分辨率为250m、500m

和1000m，扫描宽度为2330km。在对地观测过程中，每秒可同时获得6.1兆比特的来自大气、海洋和陆地表面信息，每日或每两日可获取一次全球观测数据。

本项目研究的数据来源于NASA公开下载的MODIS 1B数据。本研究中用1000m分辨率的数据，主要在于1000m分辨率的数据包含了36个波段全部数据，而且其分辨率大小适中。所获取的数据是2009年8月1日到2009年8月30日白天数据，范围覆盖四川省。

2. 1. 2气象数据源

气象数据主要有四川省12个降雨量观测站2000年到2009年每月3旬的气象观测数据，站点分布如下图所示：

图2-1气象数据站点分布图

2. 1. 3 其它数据源

四川省省界来源于国家基础地理信息中心（National Geomatics Center of China），另外还有中国区域的数字高程模型（DEM，分辨率为1km，来源于NASA）、投影全部采用埃尔伯斯等积圆锥投影（Albers Conical Equal Area），这些基础数据用于生成旱情监测系统的基础数据库。

2.2GDAL库的配置

GDAL库默认不支持项目所用的HDF文件格式。要获得支持HDF文件格式的GDAL，需要向GDAL源代码中，添加HDF文件驱动重新编译。配置过程如下：

a) 从官网下载开源库GDAL (1.6版本)和HDF4函数库(42r4版本)。HDF4

函数库放在”C:\GDAL1.6\HDF4”目录下。

b) 打开gdal文件夹下的nmake.opt，修改GDAL_HOME = "C:\GDAL1.6"，把路径改到需要把gdal安装的地方。

c) 添加HDF4驱动。编辑gdal根目录下的nmake.opt，找到“# Uncomment the following and update to enable NCSA HDF Release 4 support.”这一行。把下面两行前面的#去掉，然后改成：

HDF4_DIR = C:\GDAL1.6\42r4-win\

HDF4_LIB = $(HDF4_DIR)\dll\hd424m.lib $(HDF4_DIR)\dll\hm424m.lib\

.$(HDF4_DIR)\lib\hd424.lib $(HDF4_DIR)\lib\hm424.lib

完成后保存并关闭。

d) 打开控制台，输入: ”D:\Program Files\Microsoft Visual Studio .NET 2003\Vc7\bin\vcvars32.bat" 注册vc的编译环境。依次运行：

nmake /f makefile.vc

nmake /f makefile.vc install

nmake /f makefile.vc devinstall

如果编译不成功，这时有可能是GDAL代码自带的BUG。需要修改

gdal\frmts\leveller\levellerdataset.cpp文件171行：

{ “?, kPI / 180.0, UNITLABEL_DEGREE }

将“ “? ”修改为“ ”? “ ”，保存后重新编译。

e)编译好后，把"C:\GDAL1.6\bin"添加到系统环境变量中。并在项目预定义内容中添加以下内容：

#include "C:\\GDAL1.6\\include\\gdal_priv.h"

#include "C:\\GDAL1.6\\include\\gdal.h"

#include "C:\\GDAL1.6\\include\\cpl_string.h"

#include "C:\\GDAL1.6\\include\\cpl_conv.h"

#pragma comment (lib,"C:/GDAL1.6/lib/gdal_i.lib") //编译时链接GDAL的lib库

这样在项目中就能使用支持HDF格式的GDAL库了。

3数据处理原理及实现

3.1遥感数据处理

干旱是降水时空分布不均，长时间缺乏降水而形成的气象现象。这样就会造成局部地区作物生长过程中因供水不足，阻碍作物的正常生长而发生的水量供应不足的现象。

一般来讲，植被的生产状况主要与水分有关，水分供应程度便成了作物生长的关键因素，水分供应充足，植被生长良好，反之生长变差。植被指数的时空变化与土壤水分状况有一定的相关性，因此植被指数可以作为干旱监测的因素。植被指数是卫星传感器不同通道探测数据的线性或非线性组合，能够反映绿色植物生长和分布的特征指数。本项目采用十分常用的归一化植被指数

NDVI(Normalized Difference Vegetation Index)，其定义为：

NDVI=(CH1-CH2)/(CH1+CH2) (1) 其中，CH1和CH2分别表示近红外波段和红光波段的反射率。

地表温度(Land Surface Temperature, LST)也可用于干旱监测。地表温度是控制地球表面大多数物理、化学和生物过程的参数之一。影响地表温度变化的因素也比较多，比如地表湿度、气温、光照强度、地表材质（比如是草坪还是裸露土地，还是水泥地面，或者是沥青地面）等。

综合考虑上述的两个干旱监测因子，本项目采用植被供水指数VSWI 作为遥感干旱监测因子，其定义为：

VSWI=NDVI/Ts (2)

VSWI 的物理意义：当植被供水正常时，卫星遥感的植被指数在一定的生长期内保持在一定的范围，而卫星遥感的植被冠层温度也保持在一定的范围，如果遇到干旱，植被供水不足，一方面植被的生长受到影响，卫星遥感的植被指数将降低，另一方面当植被受旱胁迫时，为减少水分损失，叶面的气孔会部分关闭，导致叶面温度的升高，从而植被冠层温度将升高。因此可以用植被供水指数来评价旱情状况2。

遥感干旱监测模块流程图如下所示：

图3-1遥感干旱监测计算流程图

2. 1. 1遥感数据的预处理

为了减小数据量的大小，MODIS 把接收的数据转转换为整型数据存储。要获取真实的反射亮度值和辐射亮度值，就需要对获取的MODIS 数据进行一个转换。

2农业科学院农业资源与农业区划研究所. 《全国旱情监测改进方法》.

MODIS 第1和第2波段

植被指数NDVI 地表比辐射率ε31和ε32 大气水汽含量

MODIS 第2和19波段 大气透过率τ31(θ)和

τ32(θ)

MODIS 第31和32波段 分裂窗算法

地表温度LST

亮度温度T 31和T 32 热辐射强度I 31和I 32 植被供水指数VSWI

几何校正 遥感干旱监测指数图

用以下公式可得传感器所获得反射亮度值或辐射亮度值：

错误！未找到引用源。 (3)

式中：offsets 为偏移量，在HDF 文件中，其名称为reflectance_offsets 或radiance_offsets ；scales 为缩放比，在HDF 文件中名称为rflectance_scales 或radiance_scales 。

2. 1. 2云检测

由于遥感数据在有云的情况下，不能反映地表真实情况。因此这里需要进行云检测。白天判断是否有云条件如下：

0.650.86120.650.8612(0.9) or (265) or (0.7 and 285)T K T K ρρρρ+><+><

其中，错误！未找到引用源。0.65和错误！未找到引用源。0.86分别为MODIS 第一、

第二波段反射率值，T12为MODIS 第32波段亮度温度值。

2 .1. 3归一化植被指数（NDVI ）的计算

植被的生长状况主要与水分有关，水分的供应程度便成了作物生长的关键因素，水分供应充足，植被生长良好，反之生长变差。因此植被的时空变化与土壤水分状况有一定的相关性，因此植被指数可以用于监测对作物生长不利的水分胁迫环境，达到对干旱程度的监测。目前最常用的是归一化植被指数，其定义为：

错误！未找到引用源。 (4)

其中，错误！未找到引用源。1、错误！未找到引用源。2分别表示下垫面的红光波段和近红外波段的反射率，即MODIS 数据的第一波段和第二波段的反射率值。

2 .1. 4地表温度（LST ）的反演方法

现有的地表温度反演算法可以分为三大类：单通道算法、分裂窗算法、和多波段算法。其中单通道算法适合于只有一个热红外的数据。而多波段算法还在发展之中，目前没有一个简单可行的算法来进行地表温度反演。而分裂窗算法要使用两个相近的热红外波段遥感数据来进行温度反演。MODIS 数据中31和32波

段就符合分裂窗算法的要求。因此我们选择分裂窗算法。而在众多的分裂窗算法中，就简单和精度考虑，我们选择Qin et al.两因素地表反演模型。

Qin et al.提出的两因素模型是根据星上亮度的线性组合来反演地表温度。其公式如下：

错误！未找到引用源。(5)

T s为地表温度，T31和T32分别是MODIS第31和32波段的亮度温度（由其DN值计算得到），A0，A1和A2是分裂窗算法的参数，分别定义如下：

错误！未找到引用源。(6)

错误！未找到引用源。(7)

错误！未找到引用源。(8)

其中a31，b31，a32，b32是常量，由Qin et al算法可得，在地表温度0~50℃范围内，这些常量可以取a31=-64.60363，b31=0.440817，a32=-68.72575，

b32=0.473453上述公式的中间参数计算如下：

错误！未找到引用源。(9)

错误！未找到引用源。(10)

错误！未找到引用源。(11)

错误！未找到引用源。(12)

错误！未找到引用源。(13)

错误！未找到引用源。(14)

其中i是指MODIS的第31和32波段，分别为i=31或32；错误！未找到引用源。是视角为错误！未找到引用源。的大气透过率；错误！未找到引用源。是波段i的地表比辐射率。

根据MODIS图像的DN值可以计算地表热辐射强度，其公式如下：

错误！未找到引用源。(15)

其中，I i是MODIS第i波段的热辐射强度，N Di是第i波段的DN值，D Ri和D RSi分别是第i波段的辐射常量，在MODIS的头文件中。

根据普朗克函数可以求解星上亮度温度，计算公式如下：

(16)

式中：T i是MODIS第i波段的亮度温度，即式(3)的T31和T32；I i是MODIS 第i波段的热辐射强度，由式(13)得出，λi是第i波段的有效中心波长；其值可取λ31=11.03μm、λ32=12.02μm；C1和C2分别是第1是第2光谱常量；

大气透过率τi(θ)是计算地表温度的基本参数，通常是通过大气水分含量来估计。在MODIS数据中，我们通过第2和19波段来反演大气水分含量，然后通过大气水分含量与大气透过率之间的函数估计大气透过率，其可能的大气水分含量用下式估计：

错误！未找到引用源。(17) 式中：w是大气水分含量；ρi表示MODIS第i波段的地面反射率。

下表3-1为天顶视角为10度的大气透过率估计方程

水分含

量

大气透过率估计方程SEE R2 F g/cm2

τ31(10)=0.99513-0.08082w 0.0044 0.9914 804.4

0.4-2.0

τ32(10)=0.99377-0.11370w 0.0055 0.9932 1028.7

2.0-4.0 τ31(10)=1.08692-0.12759w 0.0025 0.9992 1155

3.0

τ32(10)=1.07900-0.15925w 0.0008 0.9999 173498.3

τ31(10)=1.07268-0.12571w 0.0026 0.9991 9921.6

4.0-6.0

τ32(10)=0.93821-0.12613w 0.0059 0.9955 1992.4

表3-1 天顶角为10°的大气都过率估计方程

大气透过率还受遥感视角和大气剖面温度的影响，故要对大气透过率进行视角和温度校正。

下表3-2为大气透过率的温度校正函数

波段温度校正函数温度区间

MODIS 31 δτ31(T)=0.08 T>318K

δτ31(T)=-0.05+0.00325(T31-278) 278

δτ31(T)=-0.05 T<278K

MODIS32 δτ32(T)=0.095 T>318K

δτ32(T)=-0.065+0.004(T32-278) 278

δτ32(T)=-0.065 T<278K

表3-2 大气透过率的温度校正函数

大气透过率的遥感视角校正函数如下：

δτ31(θ)=-0.00247+(2.3652?10-5)θ2(18)

δτ32(θ)=-0.00322+(3.0967?10-5)θ2(19) 天顶视角可以用下式简单估计：

θ＝Va*|D0－D i|(20)

其中Va是MODIS卫星高度的星下像元视角，Va=0.0812706°；D0是星下像元所在的列号；D i是像元i所在的列号。

经过改正后的大气透过率估计方法如下（Qin et al 2001a）：

τ31(θ)=τ31(10)+δτ31(T)-δτ31(θ) (21)

τ32(θ)=τ32(10)+δτ32(T)-δτ32(θ) (22) 从式（11）知，要得出地表温度还需知道地表比辐射率。MODIS图像的地表比辐射率可以用下式估计（Qin et al . 2004）:

εi=P v R vεiv+(1-Pv)R sεis＋dε(23)

式中εi是MODIS图像第i(i=31,32)波段的地表比辐射率；εiv和εis分别是植被和裸土在第i波段的结表比辐射率，分别取ε31v=0.98672，ε32v=0.98990，ε31s= 0.96767，ε32s=0.97790。P v是像元的植被覆盖率，通过植被指数估计(见下面)；dε是热辐射相互作用校正，R v和R s分别是植被和裸土的辐射比率，建立它们与植被覆盖度之间的关系如下(Qin et al. 2004)：

R v=0.92762+0.07033P v (24)

R s=0.99782+0.08362P v(25) 其中P v为归一化植被指数，可由式（2）得出。

我们根据Sobrino et al. (2004)的研究提出如下经验公式来估计dε：

当P v=0或者P v=1时，dε最小，为dε＝0.0

当0

当1>P v>0.5时，dε＝0.003796(1-Pv)

当P v=0.5时，dε最大，为dε＝0.001898 (26) 至此，由式（5）到式（26）可以算出地表温度。

2 .1. 5植被供水指数的计算

利用计算得出的植被指数NDVI和地表温度LST可以得出植被供水指数，公式如下：

VSWI=NDVI/Ts (27)

式中：Ts是地表温度。为了能和气象观测数据综合，这里需要把植被供水指数标准化，公式如下：

错误！未找到引用源。(28)

式中：SDI是标准化植被供水指数，范围为0-100，VSWI d是最干旱旱的VSWI，VSWI w是最湿润的VSWI。

2 .1. 6几何校正

几何校正的处理流程一般包括4 个步骤:

1、控制点的选取和投影变换及输出范围的确定;

2、求输出图像空间到输入图像空间的逆变换函数;

3、逐个像素进行几何位置变换;

4、像素灰度值内插计算。

2 .1. 7利用GDAL读取MODIS数据的方法

利用软件编程读取MODIS数据是实现上述计算的前提，本项目采用开源库GDAL完成对MODIS数据的提取。这里介绍GDAL库读取MODIS数据的方法。详细代码见附录。

GDAL所支持的一些栅格数据格式如IMG、TIFF，在GDAL的抽象数据模型中，每个文件仅有一个数据集。而一个通过HDF文件存储的MODIS数据却包含多个数据集。这些数据集放在抽象模型中的子数据域，而属性信息则存放在相应子数据在元数据中。因此是利用库GDAL提取HDF文件中的波段数据可以分为以下三步。

第一步打开HDF文件和指定的子数据集。

HDF文件含有多个子数据集，因此要打开指定的子数据集，可以分为两步：第一步，打开文件，获取子数据名列表。第二步，根据第一步中获得的列表，选择并打开指定的数据集。代码如下：

GDALAllRegister(); //注册已有的文件格式驱动

const char* filename=”C:/test.hdf”;

GDALDataset* fileDataset=(GDALDataset*)GDALOpen(filename, GA_ReadOnly)//以只读方式打开文件

char **subList=GDALGetMetadata((GDALDatasetH)file Dataset,”SUBDATASETS”);//获取子数据列表

此时获取的子数据集列表中每个列表项以下面形式表示：

SUBDATASET_n_NAME=HDF4_SDS:subdataset_type:file_name:subdataset_in dex

其中subdataset_type是预定义的HDF文件名，file_name是输入的文件名，subdataset_index是subdataset序号。

访问子数据集时，子数据集的访问名字自能是“=”后面的内容，因此这里需要一个处理去掉前面的部分。代码如下：

string subDataset1 = subList[0]; //以第一个数据集为例

subDataset1Name =

subDataset1Name.substr(subDataset1Name.find_first_of(“=”)+1);

GDALDataset*subDataset=(GDALDataset*)GDALOpen(subDataset1Name.c_str(),GA_Rea dOnly); //获取子数据集指针

第二步获取子数据集属性信息。

数据集属性信息主要有数据集名字，数据大小，数据集中波段总数，波段数据位移量，波段数据缩放量等。在GDAL的抽象数据模型中，这些属性信息存放在子数据集得元数据项里。下面以读取数据辐射位移量（radiance_scales）为例说明。

Const char*

metadata=GDALGetMetadataItem((GDALDatasetH)subDataset,“radiance_scales”,NULL );//获取数据集元数据项”radiance_scales”的信息

如果上述所获得信息是与波段有关的属性，而一般一个数据集具有多个波段，则上面所得的结果是所有波段的相应信息的集合，并且每个波段的属性之间以‘，’隔开，这时，为以后使用，还需分离上述所得信息。代码如下：

string metadataStr=metadata;

vector attributeDataStr;//储存分离后的波段属性

string tmpDataStr;//用于储存临时信息

for (string::size_type i=0;i

if (metadataS tr[i]==’,’){

attributeDataStr.push_back(tmpDataStr); //储存

tmpDataStr=””;

continue;}

tmpDataStr=tmpDataStr+metadataStr[i];}//遍历全部字符，寻找’,’并删除，并分割字符串

第三步读取波段数据。

在获取子数据集得基础上，可以通过下面代码获取该子数据集里波段数据指针。

GDALRasterBand *poBand=subDataset ->GetRasterBand(bandIndex); //bandIndex是波段顺序。

GDAL库提供多个读取波段数据的函数，最常用的函数是GDALRasterBand::RasterIO()，这个函数能根据情况自动的执行数据转换，以及对数据窗口放大和缩小。对于波段数据，这个函数横向读取数据的效率远高于纵向读取的效率。因此本文读取波段数据时采用以行为单位读取。具体代码如下：

int nXSize=poBand->GetXSize();

int nYSize=poBand->GetYSize();//获取波段数据的长宽

double *pafScanline; //定义行缓冲区

vector bandData; //定义储存波段数据容器

for (int i=0;i

pafScanline = (int *) CPLMalloc(sizeof(double)*nXSize);

bandData.push_back(pafScanline);

pafScanline=NULL;} //申请存放容器

for (int j=0;j

poBand->RasterIO( GF_Read, 0,j, nXSize, 1, bandData[j]

, nXSize, 1, GDT_Float32, 0, 0 );} //以行为单位读取数据

通过以上步骤，可以获得MODIS数据中特定的波段数据，以及其相关的属性数据。完成读取工作后，可以通过GDALClose()函数关闭打开的子数据集和HDF文件。

3.2气象数据处理与降水空间插值法比较

3.2.1 气象数据的预处理