打开文件
在打开 GDAL 支持的光栅数据存储之前,需要注册驱动程序。每种支持的格式都有一个驱动程序。
通常情况下,可以通过 GDALAllRegister() 函数来完成此操作,该函数尝试注册所有已知的驱动程序,
包括使用 GDALDriverManager::AutoLoadDrivers() 自动加载的那些 .so 文件。如果对于某些应用程序,
有必要限制驱动程序集,则可以查看 gdalallregister.cpp 中的代码。在导入 gdal 模块时,
Python 会自动调用 GDALAllRegister()。
一旦驱动程序已经注册,应用程序应该调用独立的 GDALOpen() 函数来打开数据集,
传递数据集的名称和所需的访问权限(GA_ReadOnly 或 GA_Update)。
#include "gdal_priv.h"
#include <errno.h>
int main(int argc, const char* argv[])
{
if (argc != 2) {
return EINVAL;
}
const char* pszFilename = argv[1];
GDALDatasetUniquePtr poDataset;
GDALAllRegister();
const GDALAccess eAccess = GA_ReadOnly;
poDataset = GDALDatasetUniquePtr(GDALDataset::FromHandle(GDALOpen( pszFilename, eAccess )));
if( !poDataset )
{
...; // handle error
}
return 0;
}
此代码使用C++编写,包括用于处理栅格数据和错误处理的gdal_priv.h头文件和errno.h头文件。 主函数接受两个参数:传递给程序的参数数量(argc)和包含这些参数的字符串数组(argv)。 代码检查argc是否不等于2,这意味着程序没有使用正确数量的参数调用。在这种情况下,函数返回错误代码EINVAL,表示无效的参数。 如果argc等于2,则代码从argv数组中检索第二个参数,该参数应包含栅格数据集的文件名。 代码然后创建GDALDatasetUniquePtr对象,这是用于管理GDAL数据集的智能指针。调用GDALAllRegister()函数注册所有GDAL数据集驱动程序。 接下来,代码通过创建具有GA_ReadOnly标志的常量GDALAccess对象来设置数据集的期望访问级别为“只读”。 然后,代码尝试使用GDALOpen函数打开栅格数据集,将文件名和访问级别作为参数传递。如果无法打开数据集,则代码适当处理错误(未在代码片段中显示),并返回0表示成功。 总的来说,此代码从给定的文件名初始化并打开一个只读GDAL数据集,并在提供无效的命令行参数数量时返回错误代码。
请注意,如果GDALOpen()返回NULL,这意味着打开失败,并且错误消息已经通过CPLError()发出。
如果您想控制错误如何报告给用户,请查看CPLError()文档。一般来说,GDAL使用CPLError()进行错误报告。
另外,请注意,pszFilename不一定是物理文件的名称(尽管通常是这样)。它的解释是与驱动程序相关的,可能是URL,
带有在末尾添加的附加参数的文件名,控制打开或几乎任何内容。请尽量不要将GDAL文件选择对话框限制为仅选择物理文件。
获取数据集信息
如Raster Data Model中所述,GDALDataset包含一组光栅波段,所有波段都与同一区域相关,并具有相同的分辨率。
它还具有元数据、坐标系统、地理参考转换、光栅大小和各种其他信息。
在特定但常见的“北朝向”图像且没有旋转或剪切的情况下,地理参考转换Geotransform Tutorial 如下:
adfGeoTransform [0] /*左上角 x */
adfGeoTransform [1] /* w-e(东-西)像素分辨率 */
adfGeoTransform [2] /*0*/
adfGeoTransform [3] /*左上角 y */
adfGeoTransform [4] /*0*/
adfGeoTransform [5] /* n-s(南-北)像素分辨率(负值)*/
在一般情况下,这是一个仿射变换。
如果我们想打印一些有关数据集的常规信息,我们可以执行以下操作:
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] );
}
该代码用于从GDAL(地理空间数据抽象库)栅格数据文件中提取元数据。
第一行声明了一个包含六个双精度值的数组。
接下来的三行将打印出驱动程序描述、栅格大小和投影(如果有的话)。
if语句检查栅格数据是否有有效的地理转换,如果有,就使用adfGeoTransform数组中的值打印出原点和像素大小。
总的来说,该代码提供了关于栅格数据文件的基本信息,如其大小、投影和地理转换值。
获取栅格波段
目前,通过GDAL访问栅格数据是逐个波段进行的。此外,有元数据、块大小、颜色表以及其他各种基于波段的信息可用。
下面的代码从数据集中获取一个GDALRasterBand对象(从1到GDALRasterBand::GetRasterCount()编号),
并显示一些相关信息。
GDALRasterBand *poBand;
int nBlockXSize, nBlockYSize;
int bGotMin, bGotMax;
double adfMinMax[2];
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() );
这是一个利用C++中的GDAL库提取栅格数据集信息的程序。
首先声明一个指向GDALRasterBand对象的指针,该对象将表示栅格数据集的单个波段。
接下来,声明块大小、最小值和最大值的变量,以及用于存储最小值和最大值的数组。
接下来的几行代码提取有关GDALRasterBand对象的信息,例如其块大小、数据类型和颜色解释。
然后,使用printf将此信息打印到控制台。
然后,使用poBand->GetMinimum()和poBand->GetMaximum()提取该波段的最小值和最大值。
然而,这些值可能未设置,因此使用bGotMin和bGotMax变量检查它们是否成功提取。如果没有,
则调用GDALComputeRasterMinMax()函数计算最小值和最大值。
然后,使用printf将最小值和最大值打印到控制台。
最后,该代码检查是否有该波段的任何概览,使用poBand->GetOverviewCount()检查并打印概览的数量。
如果该波段有颜色表,则使用poBand->GetColorTable()检查并打印表中颜色条目的数量。
阅读栅格数据
有几种读取栅格数据的方法, 但最常见的是通过GDALRasterBand :: RasterIO()方法。
此方法将自动处理数据类型转换, 上 /下采样和窗口处理。以下代码将读取第一行扫描数据到一个相似大小的缓冲区,
并在操作过程中将其转换为浮点数。
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 );
这段代码使用GDAL库从光栅数据集中读取一个浮点数值的单行。
首先,声明一个浮点数数组指针pafScanline。然后,使用GetXSize()方法获取波段的大小(列数),
并将其存储在名为nXSize的整数变量中。使用CPLMalloc()函数为浮点数数组分配内存,其大小计算为列数和浮点数大小的乘积。
接下来,调用RasterIO()函数从波段中读取值。它需要几个参数:
- GDALAccess枚举值GF_Read表示正在读取光栅数据集。
- 要读取的区域定义为从(0,0)开始,延伸到nXSize列和1行的单行。
- 读取值的目标缓冲区为pafScanline,列数和行数分别指定为nXSize和1。
- 光栅数据集中值的数据类型为GDT_Float32,表示32位浮点数值。
- 最后两个参数设置为0。
总之,这段代码将从光栅数据集中读取一个浮点数值的单行到一个浮点数数组中。
当不再使用pafScanline缓冲区时,应该使用CPLFree()释放它。
RasterIO调用需要以下参数。