【C/C++】超大遥感影像读取和存储 GDAL

说明：本文章转载自：http://zhan.renren.com/chinalee?gid=3602888497997597705&checked=true#nogo

GDAL栅格图像操作

    GDAL是一个操作各种栅格和矢量（由ogr这个库实现）地理数据格式的开源库。包括读取、写入、转换、处理各种栅格和矢量数据格式（有些特定的格式对一些操作如写入等不支持）。

即使不是进行地理遥感方面的应用研究，GDAL也是一个非常有用的库，因为它可以支持大量我们常见的图像数据，比如jpg，gif之类的。完整的格式清单可以到此链接查看http://www.gdal.org/formats_list.html。而且已经有包括GoogleEarth在内的很多软件都是在使用GDAL作为后台库。

本文就以VC为开发平台介绍GDAL对栅格数据的操作方法。      

Include目录是开发中需要的头文件，lib中是所需要的lib文件，在VC8中应当将其存放目录添加到目录列表中，选择菜单的“工具-选项-项目和解决方案-VC++目录”，分别在“包含文件”和“库文件”中将此两个目录添加进去。

在项目的属性页中，选择“配置属性-链接器-输入”，在“附加依赖项”中添加gdal_i-vc8.lib和gdal_id-vc8.lib两个使用GDAL中需要的静态库文件，或者在程序中添加以下两行代码也可以。

#pragma comment(lib, "gdal_i-vc8.lib")

#pragma comment(lib, "gdal_id-vc8.lib")

       Bin目录下的动态链接库文件应当放置于程序能够访问的位置，比如windows\system32中。

       此外，在程序中需要引入的头文件是gdal_priv.h。

 

现在开始用C++来对图像文件进行操作。

在打开文件之前需要首先注册所需要的驱动程序，一般来说我们可以默认注册所有支持的格式驱动，所使用的函数是GDALAllRegister()。然后就是打开文件操作。

       这里要说一个数据集的概念，也就是所谓的Dataset。在GDAL中可以说数据的核心就是Dataset，简单来说可以将Dataset就理解为图像文件，比如说一个jpeg格式的文件就是一个数据集，当然其他一些文件格式可能在一个数据集中包含多于一个文件，比如可能除了图像数据文件外还可能会有一些附加信息文件等。

在数据集下最重要组成部分就是所谓的波段band，波段可多可少，比如一个RGB真彩色的图像就有3个波段，分别代表红色绿色和蓝色波段，如果是灰度图，那可能就只有一个波段，而很多遥感图像可能就会多于3个波段。

除了波段外，数据集中还含有图像相关的坐标系投影信息，元数据信息等数据。

文件的打开使用的是GDALOpen
( const char *  pszFilename, GDALAccess  eAccess 
)，pszFilename是文件路径，eAccess 是访问权限，可以是GA_ReadOnly只读，也可以是GA_Update来对文件进行修改。比如我们以只读模式打开一个tif文件：

GDALDataset *poDataset; //数据集对象指针

GDALAllRegister();//注册驱动

     poDataset = (GDALDataset *) GDALOpen( "c:\\terra335h_EV_250_Aggr500_RefSB_b0.tif",
GA_ReadOnly );

     if(
poDataset != NULL /*检查是否正常打开文件*/)

     {

         //do
something

}

       delete poDataset; //释放资源

 

在确认poDataset不是NULL的情况下就可以对图像数据集进行操作了。

       首先来看一下有关这个图像的基本信息，比如长宽和波段数。

       cout<<"RasterXSize:"<<poDataset->GetRasterXSize()<<endl;//x方向长度

     cout<<"RasterYSize:"<<poDataset->GetRasterYSize()<<endl;//y方向长度

     cout<<"RasterCount:"<<poDataset->GetRasterCount()<<endl;//波段数量

       我的这个文件输出结果如下：

RasterXSize:5684

RasterYSize:3655

RasterCount:1

       这是个5684×3655的灰度图。

      

有关的投影信息可以由GetProjectionRef函数获得。

cout<<poDataset->GetProjectionRef()<<endl;

       输出一个字符串，是WKT格式的坐标投影信息，形如：

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235630016,AUT

HORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["deg

ree",0.0174532925199433],AUTHORITY["EPSG","4326"]]

      

有关仿射变换和GCP的信息可以通过GetGeoTransform，GetGCPCount，GetGCPProjection，GetGCPs获得。比如取得仿射变换的信息如下：

double adfGeoTransform[6];//仿射参数

poDataset->GetGeoTransform(adfGeoTransform);//获取参数

然后就可以使用以下公式来计算仿射坐标了：

Xgeo = adfGeoTransform[0] + Xpixel* adfGeoTransform[1] + Yline* adfGeoTransform[2]

Ygeo = adfGeoTransform[3] + Xpixel* adfGeoTransform[4] + Yline* adfGeoTransform[5]

       这里的坐标系是以左上角为(0,0)点，向右向下为正向。

       Xgeo表示转换后的地理x坐标，Ygeo是其相应y坐标，Xpixel是图像中像素的列坐标，Yline则是行坐标。而且Xpixel和Yline的(0,0)坐标是图像左上角像素的左上角，因此左上角像素的中心坐标是Xpixel=0.5，Yline=0.5。

 

       接下来，我们就开始进入到波段处理。波段的获取使用GetRasterBand函数，参数是波段序号，从1开始。比如还是上面那个图，只有一个波段，那要得到这个波段的操作如下：

GDALRasterBand *poBand;

poBand = poDataset->GetRasterBand( 1 );

关于波段的信息也很多，我们捡几个主要的看看。

首先也是尺寸，在一个Dataset里所有波段的尺寸都一样，也就是上面用Dataset获取的数值，就不多数了。

 

然后是数据类型。不同的图像格式的存储的数据类型是不一样的，比如我的这个tif用的是UInt16，也就是无符号的短整型，其他的用的可能就是byte型或其他类型的了。获取波段数据的数据类型信息可以如下进行：

cout<<"Data
Type:"<<poBand->GetRasterDataType()<<endl;

这里返回的信息是一个数据类型的编号，对应枚举类型GDALDataType中的一个值，在此tif的数值是2，对应GDT_UInt16，如果希望返回可读性比较强的信息，可以如下

cout<<"Data
Type:"<<GDALGetDataTypeName(poBand->GetRasterDataType())<<endl;

返回的是数据类型名的字符串：Data
Type:UInt16

 

关于波段色彩类型的情况类似数据类型，也是可以由GetColorInterpretation获取关于类型分类的GDALColorInterp枚举值，同时由GDALGetColorInterpretationName来返回名称字符串：

cout<<"Color
interpretation :"<<poBand->GetColorInterpretation()<<endl;

cout<<"Color
interpretation:"<<GDALGetColorInterpretationName(poBand->GetColorInterpretation())<<endl;

输出：

Color interpretation :1

Color interpretation:Gray

表明这个波段是个灰度图波段。如果是真彩色图可能得到的结果就是Red，Blue或Green了。

 

GetMaximum和GetMinimum分别返回的是波段数据值可能的最大值和最小值，比如UInt16类型的波段最大值是65535，最小值就是0，知道这两个值就可以将图像转换成0～255的数值用于图像的显示操作了。如：

cout<<"Max
:"<<poBand->GetMaximum()<<endl;

cout<<"Min
:"<<poBand->GetMinimum()<<endl;

 

终于要说到具体数据的读写了。关于波段数据的续写核心函数就是RasterIO。这个函数可以将图像的某一个子块读入或写入。

CPLErr GDALRasterBand::RasterIO ( GDALRWFlag  eRWFlag,
int  nXOff,
int  nYOff,
int  nXSize,
int nYSize,
void *  pData,
int  nBufXSize,
int  nBufYSize, GDALDataType  eBufType,
int  nPixelSpace,
int  nLineSpace 
)

eRWFlag：读写标志位，GF_Read读，GF_Write写。

nXOff，nYOff：读写数据块的左上角坐标

nXSize，nYSize：读写数据块的xy方向像素尺寸

pData：读写数据块缓冲区指针

nBufXSize，nBufYSize：缓冲区的像素宽度和高度

eBufType：数据类型，如GDT_UInt16等

nPixelSpace：pData中一个像素的字节大小，值为0的时候就默认是sizeof(eBufType)

nLineSpace ：pData中一行数据的字节大小，值为0的时候就默认是sizeof(eBufType)*nBufXSize

       比如说我们要读取第一行数据：

short *pafScanline;

     int  
nXSize = poBand->GetXSize();

     pafScanline = ( short *)
CPLMalloc(sizeof( short)*nXSize); //分配缓冲区空间

     poBand->RasterIO( GF_Read, 0, 0, nXSize, 1,

     pafScanline, nXSize, 1, GDT_UInt16, 0, 0 );

     CPLFree(pafScanline); //释放缓冲区空间

 

获取波段的像素数据后就可以使用GDI或GL的像素操作函数绘制图像了。另外，如果希望一次性读出整个波段的所有数据，根据GDAL文档的说法使用ReadBlock函数效率会更好一些。

至于要对有关数据和信息进行写入，一般都有读取方法对应的Set*或Write*函数可以使用，在此不再赘述。

总而言之，GDAL很好很强大。