卫星SAR数据读取与保存方法研究与软件实现

刘元廷，赵朝方，王志雄，马佑军

（中国海洋大学 海洋遥感研究所，山东 青岛 266100）

合成孔径雷达（Synthetic Aperture Radar，SAR）作为主动微波传感器可以全天时、全天候的工作，具有良好的空间分辨率，因此在环境监测、军事、气象、海洋、水文、地质、冰川等很多方面有着广阔的应用[1]，并被认为是最有效、最有潜力的卫星传感器之一。随着多种卫星SAR传感器的成功发射，SAR应用领域会更加广阔。

但是，不同SAR传感器产品的数据格式不同，造成了对SAR数据的应用非常不便。GDAL（Geospatial Data Abstraction Library）[2]是一个在X/MIT许可协议下的开源栅格空间数据转换库，几乎能够读取目前所有格式的遥感数据，很多著名的GIS类产品都使用了GDAL库，包括ESRI的ArcGIS 9.2、Google Earth和跨平台的GRASSGIS系统。因此在数据读取方面GDAL是比较理想的选择，但是它对属性读取的支持有限，所以单纯用GDAL并不能满足实际要求，还需要其他辅助库。

1 数据的读取和显示

1.1 数据读取

SAR数据读取主要采用开源GDAL库实现。默认编译的GDAL 能 够 读 取 ENVISAT、ERS、Radarsat、TerraSAR-X 的 数据格式产品，但是并不支持以HDF5格式保存的CosmoSkyMed产品，如果要支持HDF5数据格式，必须将HDF5库编译到GDAL中[3]。

在读取SAR数据前，必须打开数据并获取数据的数据集。不同的文件类型，数据集的获取方法稍有差别。

对于非HDF5格式数据集的获取方法：

GDALAllRegister（）； //注册驱动

CPLSetConfigOption （“GDAL_FILENAME_IS_UTF8”，“NO”）； //支持中文路径

GDALDataset*poDataset= （GDALDataset*）GDALOpen（lpszPathName， GA_ReadOnly）； //获取数据集

poDataset即为获取的数据集。

HDF5格式比较复杂，一个HDF5文件中可能包含一个或多个数据集。因此HDF5格式的数据集的获取方法也相对比较复杂，下面是获取所有数据集的方法[2]：

GDALRegister_HDF5（）；//注册驱动

CPLSetConfigOption（“GDAL_FILENAME_IS_UTF8”，“NO”）； //支持中文路径

GDALDataset*tempDataset；

vectordatasets； //数 据 集 类 型 的vector容器

GDALDataset*Dataset= （GDALDataset*）GDALOpen（lpszPathName，GA_ReadOnly）；

char** Ch_Sub_DS= Dataset->GetMetadata（“SUBDATASETS”）；

//获取并保存HDF5文件中所有的数据集

if（ CSLCount（Ch_Sub_DS） >0 ）

{

for（int i=0；Ch_Sub_DS[i]!=NULL；i+=2 ）

{

string tmpstr=string（Ch_Sub_DS[i]）；

tmpstr=tmpstr.substr（tmpstr.find_first_of（“=”）+1）；

tempDataset=（GDALDataset*）GDALOpen（tmpstr.c_str（），GA_ReadOnly）；

datasets.push_back（tempDataset）；

}

}

datasets为 vector容器，里面存放着获取的所有数据集。

由于1.8版本之后的GDAL默认情况下不支持中文路径，可以在编译前修改，也可以在注册完驱动后加上代码：CPLSetConfigOption（“GDAL_FILENAME_IS_UTF8”，“NO”）[4]。

GDAL在SAR数据的读取方面有着很大的优势，可以用统一的代码读取所有栅格数据，数据读取的代码有以下两种方法：

第一种：GDALReadBlock （GDALRasterBandH Band， int xoff， int yoff， void*pData）。

第 二 种 ：CPLErr GDALDataset::RasterIO （GDALRWFlag eRWFlag， int nXOff， int nYOff， int nXSize， int nYSize， void* pData， int nBufXSize， int nBufYSize， GDALDataType eBufType， int nBandCount， int * panBandMap， int nPixelSpace， int nLineSpace， eBufType* nBufXSize， int nBandSpace）；

其中eRWFlag标志是读还是写，如果是GF_Read则表示代码为读、如果是GF_Write则标志代码为写；pData是存储数据的按行存储的一维数组。两种方法各有优缺点：第一种方法速度很快，但是每次读取的行列数是固定的，由要读取数据的存储结构决定，也不能抽样读取。第二种方法速度较慢，但是读取方法灵活，不但能读取任意的行列数，而且还能按照任意比例抽样读取。在实际编程中采用了第二种方法。因为绝大多数的图像都是按行存储，因此第二种方法按行读取和存储比按列要快的多。GDAL读取数据上下颠倒，在图像显示时应该注意。

1.2 数据显示

图像显示用了另一个开源库—OpenGL[5]。OpenGL显示速度非常快而且它基本上可以显示任意格式的像素，显示的主要代码：

void glDrawPixels （GLsizei width， GLsizei height，GLenum format， GLenum type， const GLvoid*pixel）；

其中format为像素格式，主要用到的有GL_LUMINANCE（灰度格式）和GL_RGB（RGB格式）；pixel是要显示的数据的指针，也是按行存储的一维或多维数组，与GDAL读取的数据相对应。另外合理使用OpenGL的双缓存可以加快显示速度，但并不是所有的显卡都支持双缓存，尤其是Intel的集成显卡。

2 属性的读取

2.1 基本属性的读取

在使用SAR数据时，必须先获取数据的一些基本信息，如，图像的长、宽和投影信息等。在非HDF5的SAR数据格式中，GDAL可以用统一的代码读取不同数据格式的基本信息，如下所示：

GDALDataType GetRasterDataType（）； //获取数据类型

int GetRasterXSize（void ）；//获取图像长

int GetRasterYSize（void ）；//获取图像宽

int GetRasterCount（void）；//获取图像一个数据集中波段数

GDALRasterBand*GetRasterBand（ int）； //获取图像中一个数据集中的某一个波段指针

virtual const char*GetProjectionRef（void）； //获取图像的投影信息

virtual CPLErr GetGeoTransform（double*）；//获取图像的仿射变换

virtual const GDAL_GCP*GetGCPs（）；//获取图像的控制点

对于HDF5格式的数据，GDAL支持不好，除长、宽、数据类型和波段数外，其它基本信息都无法正确获取。因此必须借助其他的辅助库，具体获取方法在下节全部属性的读取中介绍。

2.2 全部属性的读取

GDAL对栅格数据属性的支持并不好，即使是已经编译进HDF5库的GDAL，也无法用统一的方法获取其基本属性，而且也不能读取HDF5的全局属性。非HDF5格式的数据也只能读取很少一部分属性。因此必须自己写代码或选择其他的库作为辅助。

2.2.1 ENVISAT和ERS全部属性的读取

ENVISAT和ERS不同SAR产品甚至相同SAR数据产品的不同模式，数据格式并不完全相同，如果自己编写代码将会比较麻烦，幸好属性的读取在其官网[6]都有源代码，下载下来稍作修改就可以编译成库或者直接添加到软件工程中使用。

2.2.2 Radarsat和TerraSAR-X全部属性的读取

Radarsat和TerraSAR-X都采用了XML格式，即数据放在Tiff数据中，属性及Tiff路径等重要信息放到XML中。两种SAR数据所采用的XML与GDAL的XML域规范并不完全相同，为了能够读取XML中的属性信息，选用开源的TinyXML库，使用递归函数实现全部属性的读取。

2.2.3 CosmoSkyMed HDF5全部属性的读取

由于HDF5的自我描述性使其结构异常复杂，导致GDAL对其支持不好，不仅全局属性无法读取，基本属性也读不全。为了能够读全其属性，需要单独编译HDF5的库，读取属性时用HDF5，读取数据时用GDAL。需要注意的是GDAL也编译进了HDF5的库，而HDF5库有些版本并不能向低版本兼容，因此，GDAL所用的HDF5库和单独编译的HDF5库版本最好相同。下面为以UTM投影为例，用HDF5库读取CosmoSkyMed HDF5格式投影信息的代码：

char string_out[80]={0}；

hid_t attr， file；

hid_t atrr_type；

herr_t ret；

file=H5Fopen （lpszPathName， H5F_ACC_RDONLY，H5P_DEFAULT）； //打开 hdf5 文件

int nzone；

CString strProject；

//获取投影方式

attr=H5Aopen_name（file，"Projection ID"）；

if（attr>=0）

{

atrr_type=H5Tcopy（H5T_C_S1）；

H5Tset_size（atrr_type， 50）；

ret=H5Aread（attr， atrr_type， string_out）；

strProject.Format（“%s”，string_out）；

ret=H5Tclose（atrr_type）；

ret=H5Aclose（attr）；

}

if（strProject== “UTM”）

{

//获取投影带号

attr=H5Aopen_name（file，"Map Projection Zone"）；

if（attr>=0）

{

ret=H5Aread（attr， H5T_NATIVE_INT， &nzone）；

ret=H5Aclose（attr）；

}

}

else//其他投影方式，在此略

{

…….

}

3 数据及属性的保存

3.1 数据保存

数据保存的格式可以是GDAL支持的任意一种。与SAR数据的读取一样，保存前也需要获取要保存的数据集，可通过下面3行代码获取：

const char*pszFormat= “GTiff”；

GDALDriver*poDriver=GetGDALDriverManager （）->GetDriverByName（pszFormat）；

GDALDataset *poDstDS = poDriver ->Create （pszDstFilename， width， height， n， eDT， NULL）；

其中pszFormat为要保存的格式；pszDstFilename为输出路径。为了后续操作的简单，保存的格式应该统一，考虑到格式的通用性，将数据保存成GeoTiff格式[7]。

对于数据的保存，GDAL也有统一的代码，和读取相同，保存也有两种方法：第一种方法为GDALWriteBlock，参数和GDALReadBlock完全相同；第二中方法和数据读取的第二种方法相同，只需要将参数eRWFlag改为GF_Write。

3.2 属性保存

对于每种数据格式属性数据都能保存到TXT中。但是由于有些属性在后续处理中可能会用到，因此最好能保存以方便及时调用。GDAL可以方便的把属性保存到XML域中，并容易读出。

3.2.1 基本属性保存

在基本属性中，投影信息、仿射系数和控制点，不仅可以读取也可以保存，下面是保存分别保存的代码：

virtual CPLErr SetProjection （const char* ）；//保存投影信息

virtual CPLErr SetGeoTransform （double* ）；//保存仿射系数

virtual CPLErr SetGCPs （ int nGCPCount， const GDAL_GCP *pasGCPList， const char*pszGCPProjection ）； //控制点

接着2.2.3节中的例子，将从HDF5数据中读出的投影信息保存，示例代码：

OGRSpatialReference oSRS；

oSRS.SetWellKnownGeogCS（"WGS84"）；

oSRS.SetUTM（nzone）；

char*chproj=NULL；//投影信息

oSRS.exportToPrettyWkt（&chproj）；

poDataset->SetProjection（chproj）；

经过这部操作后，只需要调用统一代码virtual const char*GetProjectionRef（void），就可以获取图像的投影信息。仿射系数和控制点也可以采用类似的操作。

3.2.2 其余属性保存

属性的保存：

SetMetadataItem （const char*pszName，const char*pszValue， const char*pszDomain= “” ）；

其中pszName为属性的名字，pszValue为属性的值，pszDomain一般采用默认值即可。

属性保存后的读取：

virtual const char*GetMetadataItem （const char*pszName， const char*pszDomain= “”）；

pszName为属性的名字，和保存时采用的名字必须一样，返回值即为属性的值。

4 软件实现

软件采用MFC结合GDAL等开源库实现了主要卫星SAR数据的读取、显示、放大、缩小、1:1显示、全图显示、区域选择、漫游以及移动地理坐标、图像坐标和像素值的显示，以及图像和属性的保存。另外为了更方便的操作，软件添加了鹰眼视图、直方图视图和已打开文件列表视图。以CosmoSkyMed hdf5格式数据为例，SAR图像数据显示如图1所示，全部属性保存如图2所示。

图1 SAR数据显示Fig.1 SAR data display

图2 SAR属性保存示例Fig.2 SAR attributes save example

5 结束语

针对目前主要的卫星SAR传感器，文中研究了基于GDAL和HDF5对数据的读取和保存；基于OpenGL数据显示；以及基于GDAL、HDF5和TinyXML对属性的读取和保存，并着重介绍了GDAL和HDF5联合读取HDF5文件的方法。为了后续处理的简便，最后将所有类型的数据转化为Geotiff格式，并将属性信息保存到与Geotiff文件名相同的XML文件中。软件采用面向对象结构，各功能用类实现，很容易在此基础上添加其后续功能。对自主开发SAR处理软件有一定的参考价值。

[1]Seelye Marti．海洋遥感导论[M]．蒋兴伟，等，译．北京：海洋出版社，2008：306-326．

[2]Norman Frank Warmerdam.GDAL（Geospatial DataAbstraction Library）.[EB/OL].[2012-04-20].http://www.gdal.org/.

[3]王继成，蒋狄微，谢智剑．基于GDAL的HDF文件格式栅格数据提取的研究 [J]．计算机技术与信息发展，2011（8）：62-63．WANG Ji-cheng，JIANG Di-wei，XIE Zhi-jian.Research on GDAL-based HDF raster format data extraction [J].Computing Technology and Information Development，2011，（8）:62-63.

[4]李民录．关于GDAL180中文路径不能打开的问题分析与解决[EB/OL].（2011-07-16）[2012-05-12].http://blog.csdn.net/liminlu0314/article/details/6610069.

[5]Dave Shreiner， The Khronos OpenGL ARB Working Group．OpenGL编程指南[M]．7版．李军，徐波等，译．北京：机械工业出版社，2010．

[6]European Space Agency.ESA Observing the earth[EB/OL].[2012-05-22].http://www.esa.int/esaEO/index.html.

[7]牛芩涛，盛业华．GeoTIF图像文件的数据存储格式及读写[J]．四川测绘，2004，27（3）：105-108．NIU Qin-tao，SHENG Ye-hua.The storage and read/write of GeoTIFF image file[J].Surveying and Mapping of Sichuan，2004，27（3）:105-108.