HeadFirst GDAL 3 数据读写
数据读写
本章我们就开始按顺序详细介绍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");//复制