数据读写
本章我们就开始按顺序详细介绍GDAL读写数据的过程。最后将提供完整的读写流程代码。完整代码部分中有整个创建、两影像相加、写入的流程，如果已经大致了解GDAL的读写流程，可以直接参照。
GdalDriver
首先，GDAL对每种格式提供了一个驱动 GdalDriver ， GdalDriver将对对应格式的数据进行管理，例如读取、创建、删除、重命名、复制、从已有数据创建新数据集等。所以，我们所有程序开头，都将添加GDALAllRegister() 函数，注册所有GDAL支持的数据驱动。
再次强调，所有的程序都要首先调用GDALAllRegister() 函数，否则将无法打开任何数据。
数据读取
GDAL中数据读取的步骤如下：
- 打开数据集
- 打开数据集下所需的波段
- 读取数据
下面分步讲解如何操作：
打开数据集
打开 multibands 步骤如下：使用 GDALOpen 或者 GDALOpenShared函数，传入文件名，打开数据集。 GDALOpen 和 GDALOpenShared参数一样，区别在于，在同一线程，如果是相同的文件，多个 GDALOpenShared打开的其实是同一个 GDALDataset的引用（如果在不同线程使用，为了保证线程安全，返回的将是不同的对象）。
#include "gdal_priv.h"
#include "cpl_conv.h"
//...中间的程序
//所有程序前先加上
GDALAllRegister();
//pszFilename代表文件名，GA_ReadOnly表示以只读方式打开
//也可以使用GA_Update
//GDALOpenShared和GDALOpen可以互换
GDALDataset *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
if( poDataset == NULL )
{
...;//打开失败处理
}
打开 multidatasets方式与上面稍有区别，因为有多个子数据集，所以需要获取真实数据的话，实际上是获取子数据集中的波段，可以使用gdalinfo 获取如下的 subdatasetsdomain :
Subdatasets:
SUBDATASET_1_NAME=NETCDF:"example.nc":auditTrail
SUBDATASET_1_DESC=[2x80] auditTrail (8-bit character)
SUBDATASET_2_NAME=NETCDF:"example.nc":data
SUBDATASET_2_DESC=[1x61x172] data (32-bit floating-point)
SUBDATASET_3_NAME=NETCDF:"example.nc":lat
SUBDATASET_3_DESC=[61x172] lat (32-bit floating-point)
SUBDATASET_4_NAME=NETCDF:"example.nc":lon
SUBDATASET_4_DESC=[61x172] lon (32-bit floating-point)
打开子数据集也使用 GDALOpen 或者 GDALOpenShared函数，但是传入的不是整个数据集的名称，而是需要打开的SUBDATASET_N_NAME= (_N_中N为第几个子数据集，即1 2 3 4...)后面的字符串，即子数据集名称。例如需要打开上例中data子数据集，使用如下代码：
#include "gdal_priv.h"
#include "cpl_conv.h"
//...中间的程序
//所有程序前先加上
GDALAllRegister();
//也可以使用GA_Update
//GDALOpenShared和GDALOpen可以互换
//pszFilename代表子数据集的NAME，注意双引号需要转义，GA_ReadOnly表示以只读方式打开
char * pszFilename = "NETCDF:\"example.nc\":data";
GDALDataset *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
if( poDataset == NULL )
{
...;//打开失败处理
}
打开了数据集后，两种数据接下来的处理方式都是一样的。
打开波段
GDAL中， 波段起始索引为1 ，所以循环时需要特别注意。使用GDALDataset->GetRasterBand(int BandIndex) 函数，可以打开波段。
#include "gdal_priv.h"
#include "cpl_conv.h"
//...中间的程序
//所有程序前先加上
GDALAllRegister();
GDALDataset *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );//打开文件
if( poDataset == NULL )
{
//...//打开失败处理
}
//打开数据集同上
//获取影像相关信息
int i,j,width,height,bandNum;
width = inDS->GetRasterXSize();//影像宽度
height = inDS->GetRasterYSize();//影像高度
bandNum = inDS->GetRasterCount();//影像波段数量
//将第一波段读入
double *data = new double[width*height];
GDALRasterBand *band = ds1->GetRasterBand(1);//获取第一波段，波段从1开始
//如果是在循环内，需要注意，波段从1开始；
//for(i=0;i < bandNum;i++){
// GDALRasterBand *band = ds1->GetRasterBand(i+1);//这里注意下
// ......对波段处理
//}
读取内容
使用 RasterIO函数获取波段中的具体内容，该函数的功能强大，参数较多,我们先用代码演示,然后用图来表示。
函数说明：
///CPLErr GDALRasterBand::RasterIO (
///@param eRWFlag, //读取或者写入,GF_Read或GF_Write
///@param nXOff, //起始点x坐标
///@param nYOff, //起始点y坐标
///@param nXSize, //所需读取(写入)块宽度
///@param nYSize, //所读(写)块高度
///@param * pData, //所读(写)数据,指针,
///@param nBufXSize, //一般跟nXSize一致，用于缩放图像，
/// 图像将按nBufXSize/nXsize在x尺度缩放（会自动重采样）
/// ，一般不需要调整
///@param nBufYSize, //一般跟nYSize一致
///@param eBufType, //与pData的实际类型一致,GDT_Float64
/// 代表double,其他的可以跳到定义查看
///@param nPixelSpace,
/// 设置为0为自动判断，一般设为0
/// 表示的是当前像素值和下一个像素值之间的间隔，单位是字节
/// 例:byte类型，就是1，double类型，就是8
///@param nLineSpace
/// 设置为0为自动判断，一般设为0
/// 表示的是当前行和下一行的间隔，，单位是字节
/// 例如，一行300像素，类型为int,此时就是300*4 = 1200
///@return //是否成功,成功返回CE_None ,失败返回 CE_Failure
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
)
示意图：
读取单幅影像第一波段原始数据代码示例：
#include "gdal_priv.h"
#include "cpl_conv.h"
//...中间的程序
//所有程序前先加上
GDALAllRegister();
GDALDataset *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );//打开文件
if( poDataset == NULL )
{
//...打开失败处理
}
//打开数据集同上
//获取影像相关信息
int i,j,width,height,bandNum;
width = inDS->GetRasterXSize();//影像宽度
height = inDS->GetRasterYSize();//影像高度
bandNum = inDS->GetRasterCount();//影像波段数量
//将第一波段读入
double *data = new double[width*height];
GDALRasterBand *band = ds1->GetRasterBand(1);//获取第一波段，波段从1开始
band->RasterIO(GF_Read,0,0,width,height,data,width,height,GDT_Float64,0,0);
warning
使用RasterIO函数时需要注意以下几点：
- 波段索引从一开始计算
- pData数组要有足够的大小
- pData数组数据类型要和eBufType对应
- 读取完成后必须要释放数组和关闭数据集，波段不需要关闭
关闭数据集
程序结束前，必须要关闭数据集，使用 GDALClose() 函数。
GDALClose(ds2);
数据写入
数据写入分为几个步骤
- 获取驱动
- 创建数据集
- 写入数据
- 关闭数据集
获取驱动
gdaldriver 中已经说明，不论读取还是写入， GDALDriver都是很重要的存在。如果需要输出特定格式，我们必须首先获取该格式的驱动，才能创建该格式的数据集。获取方式如下：
GDALDriver *poDriver;
//一般使用tif作为输出，如果有特殊需求，请参阅gdal文档，一般修改这里改为其他驱动名称即可，例如PNG
const char *pszFormat ="GTiff";
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);//获取特殊的驱动。
if(poDriver == NULL) {
return;
}
创建数据集
根据不同输出需要，创建数据集有两种方式：
- 输出的影像与输入的影像大小、投影都相同，只有数据不同，就可以从已经读入的数据集中创建
- 投影或大小不同，需要创建新的数据集
下面分别说明两种方式如何创建数据集，下文都是以输出GeoTiff作为示例，不同的GdalDriver 中，可以设置的参数不同，具体请参考GDAL支持格式中对每种格式的说明。
从已读入的数据集中创建
从已经读入的数据集中创建数据集，可以使用 GDALDriver 的 CreateCopy函数。
CreateCopy 函数原型如下:
//从已有的dataset中创建
//@param pszFilename, 输出文件名
//@param poSrcDS, 已有的dataset
//@param bStrict, TRUE表示严格等价，一般设置为FALSE，表示拷贝副本用作编辑
//@param papszOptions, 参数设置，具体可以参考每个driver的文档
//@param pfnProgress 回调函数指针
//@param *pProgressData 传入回调函数的数据
GDALDataset * GDALDriver::CreateCopy (
const char * pszFilename,
GDALDataset * poSrcDS,
int bStrict,
char ** papszOptions,
GDALProgressFunc pfnProgress,
void * pProgressData
)
使用 CreateCopy 函数创建数据集代码如下：
GDALDataset *OutDs;
char **papszOptions = NULL;//设置压缩、存储方式等，每种格式不同，可以直接设置为NULL
//也可以根据需要设置,例如：
//papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
//papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
//创建与ds1相同的坐标信息、相同大小的文件
OutDs = poDriver->CreateCopy(OutPut1,ds1,FALSE,papszOptions,NULL,NULL);
直接创建
直接创建数据集比从已有数据集中创建会多出创建投影的步骤：
创建数据集
直接集中创建数据集，可以使用 GDALDriver 的 Create 函数。
Create 函数原型如下:
//Create文件
//@param pszFilename, 输出文件名
//@param nXSize, 文件宽度
//@param nYSize, 文件高度
//@param nBands, 波段数
//@param eType, 数据类型
//@param papszOptions 参数设置，具体可以设置的属性参考每个driver的文档
//@return dataset
GDALDataset * GDALDriver::Create (
const char * pszFilename,
int nXSize,
int nYSize,
int nBands,
GDALDataType eType,
char ** papszOptions
)
使用 Create 函数创建数据集代码如下：
GDALDataset *OutDs;
///设置：
char **papszOptions = NULL;
papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
///...
///创建512*512，单波段的，byte类型的数据
OutDs = poDriver->Create( pszDstFilename, 512, 512, 1, GDT_Byte,
papszOptions );
创建投影
在 coordinatesystem 和 affinegeotransform中已经介绍过GDAL中数据如何投影，下面我们将用代码示例：
///设置仿射变换参数
double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
OutDs->SetGeoTransform( adfGeoTransform );
//设置投影
OGRSpatialReference oSRS;
char *pszSRS_WKT = NULL;
oSRS.SetUTM( 11, TRUE );
oSRS.SetWellKnownGeogCS( "NAD27" );
oSRS.exportToWkt( &pszSRS_WKT );
poDstDS->SetProjection( pszSRS_WKT );
CPLFree( pszSRS_WKT );//使用完后释放
写入数据
写入数据与读取数据基本相同，都是使用 GDALRasterBand 的 RasterIO函数，只是第一个参数变为 GF_Write ，下面就直接给出代码，不再另作说明。
float * outData = new float [512*512];
//对outData进行计算、处理
/************************************************************************/
/*********************************处理************************************/
/************************************************************************/
GDALRasterBand *outBand = OutDs->GetRasterBand(1);
outBand->RasterIO(GF_Write,0,0,512,512,outData,512,512,GDT_Float32,0,0);
关闭数据集
数据集必须关闭后，数据才会写入到输出文件中，否则数据将在缓存中，使用GDALClose() 函数关闭数据集。
GDALClose(OutDs);
完整代码
#include "gdal_priv.h"
int main(int argc, char* argv[])
{
GDALAllRegister();//注册类型，读取写入任何类型影像必须加入此句
//以只读方式打开文件
GDALDataset *ds1 = (GDALDataset *) GDALOpen("Input1.tif", GA_ReadOnly);
if(ds1 == NULL) {
printf("不能打开第一个文件，请检查文件是否存在！");
return -1;
}
int WidthAll = ds1->GetRasterXSize();
int HightAll = ds1->GetRasterYSize();
int BandNum = ds1->GetRasterCount();
float** InBuf1 = new float*[BandNum];
int i,j,k;
//分波段读取
for(i=0; i<BandNum; i++) {
InBuf1[i] = new float[WidthAll*HightAll];
GDALRasterBand *band = ds1->GetRasterBand(i+1);//获取波段，波段从1开始
band->RasterIO(GF_Read,
0,
0,
WidthAll,
HightAll,
InBuf1[i],
WidthAll,
HightAll,
GDT_Float32,
0,
0);
}
//打开第二个文件，假设第一个文件和第二个文件大小、波段数相同
//以只读方式打开文件
GDALDataset *ds2 = (GDALDataset *) GDALOpen("Input2.tif", GA_ReadOnly);
if(ds2 == NULL) {
printf("不能打开第二个文件，请检查文件是否存在！");
return ;
}
//第二个文件中的数据
float** InBuf2 = new float*[BandNum];
for(i=0; i<BandNum; i++) {
InBuf2[i] = new float[WidthAll*HightAll];
GDALRasterBand *band = ds2->GetRasterBand(i+1);
band->RasterIO(GF_Read,
0,
0,
WidthAll,
HightAll,
InBuf2[i],
WidthAll,
HightAll,
GDT_Float32,
0,
0);
}
//写入的数据(可以最后创建写入文件,需要先创建写入数据)
float** outBuf = new float*[BandNum];
for(i=0; i<BandNum; i++) {
outBuf[i] = new float[WidthAll*HightAll];
}
////////////////////////////////////////////////////
/////////////////////处理///////////////////////////
////////////////////////////////////////////////////
//简单相加示例，实际上，这部分最好写入函数，或者类中单独封装起来，方便使用
for(j=0; j<BandNum; j++) {
for(k=0; k<WidthAll*HightAll; k++) {
outBuf[j][k] = InBuf1[j][k]+InBuf2[j][k];
}
}
////////////////////////////////////////////////////
/////////////////////处理///////////////////////////
////////////////////////////////////////////////////
//创建写入图像：
GDALDriver *poDriver;
const char *pszFormat ="GTiff";
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
if(poDriver == NULL)
return;
GDALDataset *OutDs;
char **papszOptions = NULL;
char **papszMetadata;
//这里的参数全是指tif格式的参数，如果是其他格式，请把这里所有注释掉，或者参照文档，自行设定
//设置压缩类型，envi只认得packbits压缩.
//papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
//设置压缩比,可以不用,有些时候压缩反而更大,因为是无损的,除非图片很大
//papszOptions = CSLSetNameValue( papszOptions, "ZLEVEL", "9" );
//设置bsq或者BIP存储 bsq:BAND,bip:PIXEL
papszOptions = CSLSetNameValue( papszOptions, "INTERLEAVE", "BAND" );
OutDs = poDriver->Create( "output.tif",
WidthAll,
HightAll,
BandNum,
GDT_Float32,
papszOptions );
double geos[6];
//获取第一幅图的地理变化信息，这里如果与原图不同，请自己计算
ds1->GetGeoTransform(geos);
char *pszProjection = NULL;
char *pszPrettyWkt = NULL;
//设置为第一幅图的投影，如果与原图不同，请跳过这个if部分，到下个部分直接设置投影
//设置投影，如果需要特殊投影，请找到wkt串，自行建立投影后传入
OGRSpatialReferenceH hSRS;
if( GDALGetProjectionRef( ds1 ) != NULL ){
pszProjection = (char *) GDALGetProjectionRef( ds1 );
hSRS = OSRNewSpatialReference(NULL);
if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )
{
OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );
OutDs->SetProjection(pszPrettyWkt);
}
}
//pszPrettyWkt实际上是ogc wkt串或者是proj4
//如果指定投影的，可以在http://spatialreference.org/ 网站中搜索
//找到所需要的投影后，例如西安80 3度带，中央经线117E，就是EPSG:2384，点开后
//点击ogc wkt或者proj4，抄下来就行,proj4比较短，而且没有双引号不需要转义，比较简单
//pszPrettyWkt = "+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000
// +y_0=0 +a=6378140 +b=6356755.288157528 +units=m +no_defs";
//OutDs->SetProjection(pszPrettyWkt);
//设置坐标
OutDs->SetGeoTransform(geos);
CPLFree( pszPrettyWkt );
//写入数据到outds中
for(i=0; i<BandNum; i++) {
GDALRasterBand *outBand = OutDs->GetRasterBand(i+1);
outBand->RasterIO(GF_Write,
0,
0,
WidthAll,
HightAll,
outBuf[i],
WidthAll,
HightAll,
GDT_Float32,
0,
0);
}
//关闭dataset之后，数据才会真正写入到文件中，否则都是在缓存中！！！
GDALClose(OutDs);
GDALClose(ds1);
GDALClose(ds2);
//清理环境，删除new的数组、变量等
for(i=0; i<BandNum; i++) {
delete[] InBuf1[i];
delete[] InBuf2[i];
delete[] outBuf[i];
}
delete[] InBuf1;
delete[] InBuf2;
delete[] outBuf;
getchar();
return 0;
}
其他处理
可以使用 gdaldriver 对已有影像数据进行删除、复制、重命名等处理，可以参见GdalDriver类中的具体函数说明，下面使用一些例子来说明：
GDALAllRegister();//注册类型，读取写入任何类型影像必须加入此句
GDALDriver *poDriver;
const char *pszFormat ="GTiff";
poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
poDriver->Delete("TempFile.tif");//删除
poDirver->Rename("newName.tif","oldName.tif");//重命名
poDirver->CopyFiles("Replica.tif","original.tif");//复制
