读取栅格图像
新建栅格图像读取类
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 );
Comments NOTHING