栅格API教程

打开文件

在打开GDAL支持的栅格数据存储之前，必须注册驱动程序。每个支持的格式都有一个驱动程序。通常这是通过 GDALAllRegister() 函数，该函数尝试注册所有已知的驱动程序，包括使用 GDALDriverManager::AutoLoadDrivers() . 如果对于某些应用程序，有必要限制驱动程序集，那么检查gdalallregister.cpp中的代码可能会有帮助。导入gdal模块时，Python会自动调用GDALAllRegister（）。

一旦司机注册，申请应调用独立 GDALOpen() 函数打开数据集，传递数据集的名称和所需的访问权限（GA_ReadOnly或GA_Update）。

在C++中：

#include "gdal_priv.h"
#include "cpl_conv.h" // for CPLMalloc()
int main()
{
    GDALDataset  *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
    if( poDataset == NULL )
    {
        ...;
    }

在C中：

#include "gdal.h"
#include "cpl_conv.h" /* for CPLMalloc() */
int main()
{
    GDALDatasetH  hDataset;
    GDALAllRegister();
    hDataset = GDALOpen( pszFilename, GA_ReadOnly );
    if( hDataset == NULL )
    {
        ...;
    }

在 Python 中：

from osgeo import gdal
dataset = gdal.Open(filename, gdal.GA_ReadOnly)
if not dataset:
    ...

注意如果 GDALOpen() 返回NULL表示打开失败，并且已经通过 CPLError() . 如果要控制如何向用户报告错误，请查看 CPLError() 文档。一般来说，所有的GDAL用法 CPLError() 用于错误报告。另外，请注意pszFilename实际上不需要是物理文件的名称（尽管它通常是）。它的解释依赖于驱动程序，它可能是一个URL，一个文件名，在结尾添加了额外的参数来控制open或几乎任何东西。请不要将GDAL文件选择对话框限制为仅选择物理文件。

获取数据集信息

如中所述 栅格数据模型 ，A GDALDataset 包含栅格标注栏列表，所有标注栏都属于同一区域，并且具有相同的分辨率。它还具有元数据、坐标系、地理参照转换、栅格大小和各种其他信息。

在没有任何旋转或剪切的“北向上”图像的特殊但常见的情况下，使用地理参考变换 Geotransform教程 采用以下形式：

adfGeoTransform[0] /* top left x */
adfGeoTransform[1] /* w-e pixel resolution */
adfGeoTransform[2] /* 0 */
adfGeoTransform[3] /* top left y */
adfGeoTransform[4] /* 0 */
adfGeoTransform[5] /* n-s pixel resolution (negative value) */

在一般情况下，这是一个仿射变换。

如果要打印有关数据集的一些常规信息，可以执行以下操作：

在C++中：

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] );
}

在C中：

GDALDriverH   hDriver;
double        adfGeoTransform[6];
hDriver = GDALGetDatasetDriver( hDataset );
printf( "Driver: %s/%s\n",
        GDALGetDriverShortName( hDriver ),
        GDALGetDriverLongName( hDriver ) );
printf( "Size is %dx%dx%d\n",
        GDALGetRasterXSize( hDataset ),
        GDALGetRasterYSize( hDataset ),
        GDALGetRasterCount( hDataset ) );
if( GDALGetProjectionRef( hDataset ) != NULL )
    printf( "Projection is `%s'\n", GDALGetProjectionRef( hDataset ) );
if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )
{
    printf( "Origin = (%.6f,%.6f)\n",
            adfGeoTransform[0], adfGeoTransform[3] );
    printf( "Pixel Size = (%.6f,%.6f)\n",
            adfGeoTransform[1], adfGeoTransform[5] );
}

在 Python 中：

print("Driver: {}/{}".format(dataset.GetDriver().ShortName,
                            dataset.GetDriver().LongName))
print("Size is {} x {} x {}".format(dataset.RasterXSize,
                                    dataset.RasterYSize,
                                    dataset.RasterCount))
print("Projection is {}".format(dataset.GetProjection()))
geotransform = dataset.GetGeoTransform()
if geotransform:
    print("Origin = ({}, {})".format(geotransform[0], geotransform[3]))
    print("Pixel Size = ({}, {})".format(geotransform[1], geotransform[5]))

获取栅格标注栏

此时，通过GDAL对栅格数据的访问是一次一个波段完成的。此外，还提供了元数据、块大小、颜色表和各种其他信息，这些信息都是按波段提供的。以下代码获取 GDALRasterBand 数据集中的对象（编号1到 GDALRasterBand::GetRasterCount() )并显示一些关于它的信息。

在C++中：

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中：

GDALRasterBandH hBand;
int             nBlockXSize, nBlockYSize;
int             bGotMin, bGotMax;
double          adfMinMax[2];
hBand = GDALGetRasterBand( hDataset, 1 );
GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );
printf( "Block=%dx%d Type=%s, ColorInterp=%s\n",
        nBlockXSize, nBlockYSize,
        GDALGetDataTypeName(GDALGetRasterDataType(hBand)),
        GDALGetColorInterpretationName(
            GDALGetRasterColorInterpretation(hBand)) );
adfMinMax[0] = GDALGetRasterMinimum( hBand, &bGotMin );
adfMinMax[1] = GDALGetRasterMaximum( hBand, &bGotMax );
if( ! (bGotMin && bGotMax) )
    GDALComputeRasterMinMax( hBand, TRUE, adfMinMax );
printf( "Min=%.3fd, Max=%.3f\n", adfMinMax[0], adfMinMax[1] );
if( GDALGetOverviewCount(hBand) > 0 )
    printf( "Band has %d overviews.\n", GDALGetOverviewCount(hBand));
if( GDALGetRasterColorTable( hBand ) != NULL )
    printf( "Band has a color table with %d entries.\n",
            GDALGetColorEntryCount(
                GDALGetRasterColorTable( hBand ) ) );

在 Python 中：

band = dataset.GetRasterBand(1)
print("Band Type={}".format(gdal.GetDataTypeName(band.DataType)))

min = band.GetMinimum()
max = band.GetMaximum()
if not min or not max:
    (min,max) = band.ComputeRasterMinMax(True)
print("Min={:.3f}, Max={:.3f}".format(min,max))

if band.GetOverviewCount() > 0:
    print("Band has {} overviews".format(band.GetOverviewCount()))

if band.GetRasterColorTable():
    print("Band has a color table with {} entries".format(band.GetRasterColorTable().GetCount()))

读取栅格数据

有几种读取栅格数据的方法，但最常见的是通过 GDALRasterBand::RasterIO() 方法。此方法将自动处理数据类型转换、上/下采样和窗口。下面的代码将把第一个扫描行的数据读入一个类似大小的缓冲区，作为操作的一部分将其转换为浮点。

在C++中：

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 );

当pafScanline缓冲区不再使用时，应使用CPLFree（）释放它。

在C中：

float *pafScanline;
int   nXSize = GDALGetRasterBandXSize( hBand );
pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);
GDALRasterIO( hBand, GF_Read, 0, 0, nXSize, 1,
            pafScanline, nXSize, 1, GDT_Float32,
            0, 0 );

当pafScanline缓冲区不再使用时，应使用CPLFree（）释放它。

在 Python 中：

scanline = band.ReadRaster(xoff=0, yoff=0,
                        xsize=band.XSize, ysize=1,
                        buf_xsize=band.XSize, buf_ysize=1,
                        buf_type=gdal.GDT_Float32)

注意，返回的扫描线是string类型，包含xsize*4字节的原始二进制浮点数据。可以使用标准库中的struct模块将其转换为Python值：

import struct
tuple_of_floats = struct.unpack('f' * b2.XSize, scanline)

RasterIO调用接受以下参数。

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 )

请注意，基于eRWFlag的设置（GF-read或GF-write），相同的RasterIO（）调用用于读取或写入。nXOff、nYOff、nXSize、nYSize参数描述要读取（或写入）的磁盘上栅格数据的窗口。它不必落在平铺边界上，尽管如果这样做，访问可能会更高效。

pData是数据读入或写入的内存缓冲区。它的真正类型必须是作为eBufType传递的任何类型，例如GDT_Float32或GDT_Byte。RasterIO（）调用将负责在缓冲区的数据类型和频带的数据类型之间进行转换。请注意，当将浮点数据转换为整数RasterIO（）时，向下舍入；当将源值转换到输出的合法范围之外时，将使用最接近的合法值。这意味着，例如，读取到GDT_字节缓冲区中的16位数据将映射所有大于255到255的值，数据不会缩放！

nBufXSize和nBufYSize值描述缓冲区的大小。当以全分辨率加载数据时，这将与窗口大小相同。但是，要加载分辨率降低的概述，可以将其设置为小于磁盘上的窗口。在这种情况下，如果概览合适，RasterIO（）将利用概览来更有效地执行IO。

nPixelSpace和nLineSpace通常为零，表示应使用默认值。然而，它们可用于控制对存储器数据缓冲器的访问，例如允许读入包含其它像素交错数据的缓冲器。

关闭数据集

请记住 GDALRasterBand 对象是由它们的数据集拥有的，并且它们永远不应该被C++删除操作符破坏。 GDALDataset 可以通过呼叫关闭 GDALClose() （对于Windows用户，不建议在GDALDataset上使用delete运算符，因为在跨模块边界分配和释放内存时存在已知问题。请参见常见问题解答中的相关主题）。调用GDALClose将导致正确的清理，并刷新任何挂起的写操作。如果忘记对在更新模式下以GTiff等流行格式打开的数据集调用GDALClose，很可能会导致以后无法打开它。

创建文件的技术

如果格式驱动程序支持创建，则可以创建GDAL支持格式的新文件。使用CreateCopy（）和Create（）创建文件有两种通用技术。CreateCopy方法涉及调用格式驱动程序上的CreateCopy（）方法，并传入应复制的源数据集。Create方法涉及在驱动程序上调用Create（）方法，然后使用单独的调用显式地写入所有元数据和栅格数据。所有支持创建新文件的驱动程序都支持CreateCopy（）方法，但只有少数驱动程序支持Create（）方法。

要确定特定格式是否支持Create或CreateCopy，可以检查格式驱动程序对象上的DCAP_Create和DCAP_CreateCopy元数据。确保 GDALAllRegister() 在打电话之前 GDALDriverManager::GetDriverByName() . 在本例中，我们获取一个驱动程序，并确定它是否支持Create（）和/或CreateCopy（）。

在C++中：

#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 );

在C中：

#include "cpl_string.h"
...
const char *pszFormat = "GTiff";
GDALDriverH hDriver = GDALGetDriverByName( pszFormat );
char **papszMetadata;
if( hDriver == NULL )
    exit( 1 );
papszMetadata = GDALGetMetadata( hDriver, NULL );
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 );

在 Python 中：

fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
metadata = driver.GetMetadata()
if metadata.get(gdal.DCAP_CREATE) == "YES":
    print("Driver {} supports Create() method.".format(fileformat))

if metadata.get(gdal.DCAP_CREATECOPY) == "YES":
    print("Driver {} supports CreateCopy() method.".format(fileformat))

注意，许多驱动程序是只读的，不支持Create（）或CreateCopy（）。

使用CreateCopy（）

这个 GDALDriver::CreateCopy() 方法的使用相当简单，因为大多数信息都是从源数据集中收集的。但是，它包括用于传递特定于格式的创建选项的选项，以及用于在进行长数据集复制时向用户报告进度的选项。从一个名为pszSrcFilename的文件到一个名为pszDstFilename的新文件的简单复制，使用先前获取驱动程序的格式上的默认选项，可能如下所示：

在C++中：

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 );

在C中：

GDALDatasetH hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
GDALDatasetH hDstDS;
hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE,
                        NULL, NULL, NULL );
/* Once we're done, close properly the dataset */
if( hDstDS != NULL )
    GDALClose( hDstDS );
GDALClose(hSrcDS);

在 Python 中：

src_ds = gdal.Open(src_filename)
dst_ds = driver.CreateCopy(dst_filename, src_ds, strict=0)
# Once we're done, close properly the dataset
dst_ds = None
src_ds = None

注意，CreateCopy（）方法返回一个可写的数据集，必须正确关闭它才能完成数据集的写入并将其刷新到磁盘。在Python的情况下，当“dst_ds”超出范围时会自动发生这种情况。在CreateCopy（）调用中的目标文件名之后，bStrict选项使用的FALSE（或0）值指示，即使无法创建与输入数据集完全匹配的目标数据集，CreateCopy（）调用也应继续进行，而不会出现致命错误。这可能是因为输出格式不支持输入数据集的像素数据类型，或者是因为目标无法支持例如写入地理参考。

更复杂的情况可能涉及传递创建选项和使用预定义的进度监视器，如下所示：

在C++中：

#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 );

在C中：

#include "cpl_string.h"
...
char **papszOptions = NULL;
papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE,
                        papszOptions, GDALTermProgres, NULL );
/* Once we're done, close properly the dataset */
if( hDstDS != NULL )
    GDALClose( hDstDS );
CSLDestroy( papszOptions );

在 Python 中：

src_ds = gdal.Open(src_filename)
dst_ds = driver.CreateCopy(dst_filename, src_ds, strict=0,
                        options=["TILED=YES", "COMPRESS=PACKBITS"])
# Once we're done, close properly the dataset
dst_ds = None
src_ds = None

使用Create（）

对于不只是将现有文件导出到新文件的情况，通常需要使用 GDALDriver::Create() 方法（尽管通过使用虚拟文件或内存文件可以使用一些有趣的选项）。Create（）方法采用与CreateCopy（）类似的选项列表，但必须显式提供图像大小、带区数和带区类型。

在C++中：

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

在C中：

GDALDatasetH hDstDS;
char **papszOptions = NULL;
hDstDS = GDALCreate( hDriver, pszDstFilename, 512, 512, 1, GDT_Byte,
                    papszOptions );

在 Python 中：

dst_ds = driver.Create(dst_filename, xsize=512, ysize=512,
                    bands=1, eType=gdal.GDT_Byte)

成功创建数据集后，必须将所有适当的元数据和栅格数据写入文件。这是什么将根据使用不同，但一个简单的情况下，投影，地理转换和栅格数据在这里介绍。

在C++中：

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 );

在C中：

double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
OGRSpatialReferenceH hSRS;
char *pszSRS_WKT = NULL;
GDALRasterBandH hBand;
GByte abyRaster[512*512];
GDALSetGeoTransform( hDstDS, adfGeoTransform );
hSRS = OSRNewSpatialReference( NULL );
OSRSetUTM( hSRS, 11, TRUE );
OSRSetWellKnownGeogCS( hSRS, "NAD27" );
OSRExportToWkt( hSRS, &pszSRS_WKT );
OSRDestroySpatialReference( hSRS );
GDALSetProjection( hDstDS, pszSRS_WKT );
CPLFree( pszSRS_WKT );
hBand = GDALGetRasterBand( hDstDS, 1 );
GDALRasterIO( hBand, GF_Write, 0, 0, 512, 512,
            abyRaster, 512, 512, GDT_Byte, 0, 0 );
/* Once we're done, close properly the dataset */
GDALClose( hDstDS );

在 Python 中：

from osgeo import osr
import numpy
dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("NAD27")
dst_ds.SetProjection(srs.ExportToWkt())
raster = numpy.zeros((512, 512), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
# Once we're done, close properly the dataset
dst_ds = None