GDAL学习笔记（3）

读取栅格图像

新建栅格图像读取类

class Raster_demo {
public:
	//打开图像
	static void open(const char*, GDALDatasetUniquePtr&);
	//打印图像信息
	static void printDatasetInfo(GDALDatasetUniquePtr&);
	//抓取波段
	static GDALRasterBand* fetchRasterband(GDALDatasetUniquePtr& poDataset);
	//读取图像
	static void readRasterband(GDALRasterBand*);
};

打开图像

void Raster_demo::open(const char* pszFilename, GDALDatasetUniquePtr& poDataset) {
	GDALAllRegister();
	const GDALAccess eAccess = GA_ReadOnly;
	poDataset = GDALDatasetUniquePtr(GDALDataset::FromHandle(GDALOpen(pszFilename, eAccess)));
	if (!poDataset)
	{
		cout << "open file failed" << endl;
		return;
	}
}

获取栅格基础信息

基础信息中的仿射地理坐标变换参数很重要

void Raster_demo::printDatasetInfo(GDALDatasetUniquePtr& poDataset) {
	double        adfGeoTransform[6];
	printf("Driver: %s/%s\n",
		poDataset->GetDriver()->GetDescription(),
		poDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME));
	printf("Size is %dx%dx%d\n",
		poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
		poDataset->GetRasterCount());
	if (poDataset->GetProjectionRef() != NULL)
		printf("Projection is `%s'\n", poDataset->GetProjectionRef());
	if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None)
	{
		printf("Origin = (%.6f,%.6f)\n",
			adfGeoTransform[0], adfGeoTransform[3]);
		printf("Pixel Size = (%.6f,%.6f)\n",
			adfGeoTransform[1], adfGeoTransform[5]);
	}
}

获取波段

获取某一波段以及相关的信息。

GDALRasterBand* Raster_demo::fetchRasterband(GDALDatasetUniquePtr& poDataset) {
	GDALRasterBand* poBand;
	int nBlockXSize, nBlockYSize;
	int bGotMin, bGotMax;
	double adfMinMax[2];
	cout << "Raster Band Count = " << poDataset->GetRasterCount() << endl;
	poBand = poDataset->GetRasterBand(1);
	poBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
	printf("Block=%dx%d Type=%s, ColorInterp=%s\n",
		nBlockXSize, nBlockYSize,
		GDALGetDataTypeName(poBand->GetRasterDataType()),
		GDALGetColorInterpretationName(
			poBand->GetColorInterpretation()));
	adfMinMax[0] = poBand->GetMinimum(&bGotMin);
	adfMinMax[1] = poBand->GetMaximum(&bGotMax);
	if (!(bGotMin && bGotMax))
		GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
	printf("Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1]);
	if (poBand->GetOverviewCount() > 0)
		printf("Band has %d overviews.\n", poBand->GetOverviewCount());
	if (poBand->GetColorTable() != NULL)
		printf("Band has a color table with %d entries.\n",
			poBand->GetColorTable()->GetColorEntryCount());
	return poBand;
}

读取图像

一次按一个波段进行读取，读取前需先根据栅格数据类型和一次读取大小开辟内存

void Raster_demo::readRasterband(GDALRasterBand* poBand) {
	float* pafScanline;
	int   nXSize = poBand->GetXSize();
	pafScanline = (float*)CPLMalloc(sizeof(float) * nXSize);
	poBand->RasterIO(GF_Read, 0, 0, nXSize, 1,
		pafScanline, nXSize, 1, GDT_Float32,
		0, 0);

	//需要释放
	CPLFree(pafScanline);
}

释放资源

GDALRasterBand不能使用delete删除，会在GDALDataset释放时一并释放。使用GDALClose函数释放GDALDataset，GDALDatasetUniquePtr的删除函数就是调用的GDALClose，因此可以不管或者用reset释放。

创建栅格图像

检查相应数据格式驱动是否可用

#include "cpl_string.h"
...
    const char *pszFormat = "GTiff";
    GDALDriver *poDriver;
    char **papszMetadata;
    poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if( poDriver == NULL )
        exit( 1 );
    papszMetadata = poDriver->GetMetadata();
    if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
        printf( "Driver %s supports Create() method.\n", pszFormat );
    if( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE ) )
        printf( "Driver %s supports CreateCopy() method.\n", pszFormat );

使用CreateCopy创建

GDALDataset *poSrcDS =
(GDALDataset *) GDALOpen( pszSrcFilename, GA_ReadOnly );
GDALDataset *poDstDS;
poDstDS = poDriver->CreateCopy( pszDstFilename, poSrcDS, FALSE,
                                NULL, NULL, NULL );
/* Once we're done, close properly the dataset */
if( poDstDS != NULL )
    GDALClose( (GDALDatasetH) poDstDS );
GDALClose( (GDALDatasetH) poSrcDS );

增加一些参数设置

#include "cpl_string.h"
...
char **papszOptions = NULL;
papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
poDstDS = poDriver->CreateCopy( pszDstFilename, poSrcDS, FALSE,
                                papszOptions, GDALTermProgress, NULL );
/* Once we're done, close properly the dataset */
if( poDstDS != NULL )
    GDALClose( (GDALDatasetH) poDstDS );
CSLDestroy( papszOptions );

使用Create创建

使用create创建栅格文件（图像尺寸，波段数，数据类型必须给出）

GDALDataset *poDstDS;
char **papszOptions = NULL;
poDstDS = poDriver->Create( pszDstFilename, 512, 512, 1, GDT_Byte,
                            papszOptions );

dataset创建后，所有适当的元数据和栅格数据必须写入文件。

double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
OGRSpatialReference oSRS;
char *pszSRS_WKT = NULL;
GDALRasterBand *poBand;
GByte abyRaster[512*512];
poDstDS->SetGeoTransform( adfGeoTransform );
oSRS.SetUTM( 11, TRUE );
oSRS.SetWellKnownGeogCS( "NAD27" );
oSRS.exportToWkt( &pszSRS_WKT );
poDstDS->SetProjection( pszSRS_WKT );
CPLFree( pszSRS_WKT );
poBand = poDstDS->GetRasterBand(1);
poBand->RasterIO( GF_Write, 0, 0, 512, 512,
                abyRaster, 512, 512, GDT_Byte, 0, 0 );
/* Once we're done, close properly the dataset */
GDALClose( (GDALDatasetH) poDstDS );