gdal GDALRasterIO 具体分析

CPLErr CPL_DLL CPL_STDCALL GDALRasterIO	(	GDALRasterBandH	hBand,
		GDALRWFlag	eRWFlag,
		int	nXOff,
		int	nYOff,
		int	nXSize,
		int	nYSize,
		void *	pData,
		int	nBufXSize,
		int	nBufYSize,
		GDALDataType	eBufType,
		int	nPixelSpace,
		int	nLineSpace
	)

找到了一个虽然不是GDALRasterIO的解释方法，但是的跟它非常相似的

http://hi.baidu.com/201021170051/item/8a1ffe1c091eb9f765eabfb4

GDAL的GDALDataset和GDALRasterBand类为我们提供了功能强大的RasterIO(…)函数，它可以把图像上 指定大小的矩形象素块 以缩放的形式 按指定的数据类型 输出或输入到 用户指定大小的缓冲区中。
原型：
CPLErr GDALDataset::RasterIO(
GDALRWFlag eRWFlag, //读写标记
int nXOff, int nYOff,//相对于图像左上角顶点(从零开始)的行列偏移量 
int nXSize, int nYSize,//要读写的块在x方向的象素个数和y方向的象素列数
void * pData,//指向目标缓冲区的指针，由用户分配
int nBufXSize, int   nBufYSize,//目标块在x方向上和y方向上的大小
GDALDataType eBufType, //目标缓冲区的数据类型，原类型会自动转换为目标类型
int nBandCount, //要处理的波段数
int * panBandMap,//记录要操作的波段的索引(波段索引从1开始)的数组，若为空，//则数组中存放的是前nBandCount个波段的索引
int nPixelSpace,//x方向上两个相邻象素之间的字节偏移，默认为0,则列间的实际//字节偏移由目标数据类型eBufType确定
int nLineSpace, //y方向上相邻两行之间的字节偏移, 默认为0,则行间的实际字节//偏移为eBufType * nBufXSize
int nBandSpace   //相邻两波段之间的字节偏移，默认为0,则意味着波段是顺序结构//的，其间字节偏移为nLineSpace * nBufYSize
); 

CPLErr GDALRasterBand::RasterIO(//各参数意义同上
GDALRWFlag eRWFlag, 
int nXOff, int nYOff,
int nXSize, int nYSize, 
void * pData,//
int nBufXSize, int   nBufYSize,
GDALDataType eBufType,
int nBandCount, 
int nPixelSpace,
int nLineSpace, 
); 
上述函数的区别在于前者是对同一数据集的多个波段操作，而后者是对单个波段的操作。

基于以上的GDALRasterBand::RasterIO()函数构建影像分块的类clsDOM，并在类中做了函数SmartRead()进行缩放读取dom，处理过程为：
//打开文件获得指向该文件的GDALDataset指针
clsDOM::OpenFile();
//指定分块的大小
clsDOM::SetBlockSize();
//循环处理每一块
{
//SmartRead()内部调用RasterIO(),此时动态分配当前块的内存
clsDOM::SmartRead();
//每一块处理完之后立即释放当前块
}
//所有块处理完后，释放GDALDataset指针
clsDOM::CloseFile();
跟踪以上过程发现读取各块时内存逐次增长，虽然每次循环都释放了当前块，但是仍无法回落到处理这一块之前的值，直到最终释放GDALDataset指针内存才会回落至打开文件前的值。内存的不断增长显然对处理大文件会产生致命问题！究其原因，每执行一次RasterIO内存都会（不太慢的）增长，以致于每一块处理完后内存不彻底回落。尝试每次处理完后删除当前的RasterBand指针也无法回落。最终解决办法是处理每一块时都执行一次打开文件和关闭文件释放GDALDataset指针的操作：
//指定分块的大小
clsDOM::SetBlockSize();
//循环处理每一块
{
//处理每一块之前打开文件获得指向该文件的GDALDataset指针
clsDOM::OpenFile();
//SmartRead()内部调用RasterIO(),此时动态分配当前块的内存
SmartRead();
//每一块处理完之后立即释放当前块
//…
//当前块处理完后，释放GDALDataset指针
CloseFile();
}
这样做内存的开销极小，但频繁的打开关闭文件总是不太好。
另外，如果是数据库连接，像下面这样反复连接更是不可取的：
       GDALOpenInfo   OpenInfo(szdbConnect , GA_ReadOnly);
       GDALDataset poDataset=SDEDataset::Open(&OpenInfo);

//以下是以前的

RasterIO()可以通过指定eRWFlag参数来确定是读/写数据（GF_Read或GF_Write）。 参数nXOff/nYOff/nXSize/nYSize描述了要读的影象范围（或者是写）。同时它也可以 自动处理边界等特殊情况。

参数pData指定读/写对应的缓冲。缓冲的类型必须是eBufType中定义的， 例如GDT_Float32、GDT_Byte等。RasterIO ()会自动转换缓冲和波段的类型， 使它们一致。当数据向下转换时，或者是数据超出转换后的数据类型可以表示的范围时， 将会用最接近的数据来代替。例如一个 16位的整数被转换为GDT_Byte时，所有大于255的 值都会用255代替（数据并不会被缩放）

参数nBufXSize和nBufYSize描述了缓冲的大小。当时读写是是全部数据时， 该值和影象的大小相同。当需要对影象抽样的时候，缓冲也可以比真实的影象小。 因此，利用RasterIO()实现预览功能是很方便的。

参数nPixelSpace和nLineSpace通常被设置为0。当然，也可以使用他们来控制内存中的数据。 关闭Dataset

需要强调的一点是：GDALRasterBand对象属于相应的dataset，用户不能私自delete 任何GDALRasterBand对象。GDALDataset可以用GDALClose()关闭数据，或者是直接 delete GDALDataset对象。关闭GDALDataset的时候会进行相关的清除操作和刷新一些写操作。

后来又觉得它的整篇文档都写得很不错，所以摘录如下：

http://blog.sciencenet.cn/blog-98709-486392.html

在打开GDAL所支持的光栅数据之前需要注册驱动。这里的驱动是针对GDAL支持的所有 数据格式。通常可以通过调用 GDALAllRegister() 函数来注册所有已知的驱动，同时 也包含那些用 GDALDriverManager::AutoLoadDrivers() 从.so文件中自动装载驱动。 如果程序需要对某些驱动做限制，可以参考 gdalallregister.cpp 代码。
当驱动被注册之后，我们就可以用 GDALOpen() 函数来打开一个数据集。打开的方式 可以是 GA_ReadOnly 或者 GA_Update。

In C++:

#include "gdal_priv.h"

int main()
{
    GDALDataset  *poDataset;

    GDALAllRegister();

    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
    if( poDataset == NULL )
    {
        ...;
    }

In C:

#include "gdal.h"

int main()
{
    GDALDatasetH  hDataset;

    GDALAllRegister();

    hDataset = GDALOpen( pszFilename, GA_ReadOnly );
    if( hDataset == NULL )
    {
        ...;
    }

In Python:

    import gdal
    from gdalconst import *

    dataset = gdal.Open( filename, GA_ReadOnly )
    if dataset is None:
        ...

如果 GDALOpen() 函数返回NULL则表示打开失败，同时 CPLError() 函数产生相应的错误信息。 如果您需要对错误进行处理可以参考 CPLError() 相关文档。通常情况下，所有的 GDAL函数都通过CPLError()报告错误。另外需要注意的是pszFilename并不一定对应一个 实际的文件名（当然也可以就是一个文件名）。它的具体解释由相应的驱动程序负责。 它可能是一个URL，或者是文件名以后后面带有许多用于控制打开方式的参数。通常建议， 不要在打开文件的选择对话框中对文件的类型做太多的限制。

获取Dataset信息
如果GDAL数据模型一节所描述的，一个GDALDataset包含了光栅数据的一系列的波段信息。 同时它还包含元数据、一个坐标系统、投影类型、光栅的大小以及其他许多信息。

    adfGeoTransform[0] /* 左上角 x */
    adfGeoTransform[1] /* 东西方向一个像素对应的距离 */
    adfGeoTransform[2] /* 旋转, 0表示上面为北方 */
    adfGeoTransform[3] /* 左上角 y */
    adfGeoTransform[4] /* 旋转, 0表示上面为北方 */
    adfGeoTransform[5] /* 南北方向一个像素对应的距离 */

如果需要输出dataset的基本信息，可以这样：

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

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

In Python:

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

获取一个光栅波段
现在，我们可以通过GDAL获取光栅的一个波段。同样每个波段含有元数据、块大小、 颜色表以前其他一些信息。下面的代码从dataset获取一个GDALRasterBand对象， 并且显示它的一些信息。
In 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() );

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

In Python:

        band = dataset.GetRasterBand(1)

        print 'Band Type=',gdal.GetDataTypeName(band.DataType)

        min = band.GetMinimum()
        max = band.GetMaximum()
        if min is not None and max is not None:
            (min,max) = ComputeRasterMinMax(1)
        print 'Min=%.3f, Max=%.3f' % (min,max)

        if band.GetOverviewCount() > 0:
            print 'Band has ', band.GetOverviewCount(), ' overviews.'

        if not band.GetRasterColorTable() is None:
            print 'Band has a color table with ', \
            band.GetRasterColorTable().GetCount(), ' entries.'

读光栅数据
GDAL有几种读光栅数据的方法，但是GDALRasterBand::RasterIO()是最常用的一种。 该函数可以自动转换数据类型、采样以及裁剪。下面的代码读光栅的第1行数据， 同时转换为float保存到缓冲。
In 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 );

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

In Python:

        scanline = band.ReadRaster( 0, 0, band.XSize, 1, \
                                     band.XSize, 1, GDT_Float32 )

返回的是一个string，包含了xsize*4大小的二进制数据，是float类型指针。 可以使用python的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 )

RasterIO()可以通过指定eRWFlag参数来确定是读/写数据（GF_Read或GF_Write）。 参数nXOff/nYOff/nXSize/nYSize描述了要读的影象范围（或者是写）。同时它也可以 自动处理边界等特殊情况。

参数pData指定读/写对应的缓冲。缓冲的类型必须是eBufType中定义的， 例如GDT_Float32、GDT_Byte等。RasterIO ()会自动转换缓冲和波段的类型， 使它们一致。当数据向下转换时，或者是数据超出转换后的数据类型可以表示的范围时， 将会用最接近的数据来代替。例如一个 16位的整数被转换为GDT_Byte时，所有大于255的 值都会用255代替（数据并不会被缩放）。

参数nBufXSize和nBufYSize描述了缓冲的大小。当时读写是是全部数据时， 该值和影象的大小相同。当需要对影象抽样的时候，缓冲也可以比真实的影象小。 因此，利用RasterIO()实现预览功能是很方便的。

参数nPixelSpace和nLineSpace通常被设置为0。当然，也可以使用他们来控制内存中的数据。 关闭Dataset

需要强调的一点是：GDALRasterBand对象属于相应的dataset，用户不能私自delete 任何GDALRasterBand对象。GDALDataset可以用GDALClose()关闭数据，或者是直接 delete GDALDataset对象。关闭GDALDataset的时候会进行相关的清除操作和刷新一些写操作。

创建文件的技巧
如果相应格式的驱动支持写操作的话，则可以创建文件。GDAL有两函数可以创建文件： CreateCopy()和Create()。 CreateCopy()函数直接从参数给定的数据集复制数据。 Create()函数则需要用户明确地写入各种数据（元数据、光栅数据等）。所有支持创建 的格式驱动都支持CreateCopy()函数，但是并不一定支持Create()函数。
为了确定数据格式是否支持Create或CreateCopy，可以检查驱动对象中的DCAP_CREATE 和DCAP_CREATECOPY元数据。在使用GetDriverByName()函数之前确保GDALAllRegister() 已经被调用过。

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

In C:

#include "cpl_string.h"
...
    const char *pszFormat = "GTiff";
    GDALDriver 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 );

In Python:

    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    metadata = driver.GetMetadata()
    if metadata.has_key(gdal.DCAP_CREATE) \
       and metadata[gdal.DCAP_CREATE] == 'YES':
        print 'Driver %s supports Create() method.' % format
    if metadata.has_key(gdal.DCAP_CREATECOPY) \
       and metadata[gdal.DCAP_CREATECOPY] == 'YES':
        print 'Driver %s supports CreateCopy() method.' % format

我们可以看出有些格式不支持Create()或CreateCopy()调用。

使用CreateCopy()
GDALDriver::CreateCopy()函数使用比较简单，并且原先数据中的所有信息都被正确 的设置。函数还可以 指定某些可的选择参数，也通过一个回调函数来获得数据复制的 进展情况。下面的程序用默认的方式copy一个pszSrcFilename文件，保存 为 pszDstFilename 文件。
In C++:

    GDALDataset *poSrcDS = 
       (GDALDataset *) GDALOpen( pszSrcFilename, GA_ReadOnly );
    GDALDataset *poDstDS;

    poDstDS = poDriver->CreateCopy( pszDstFilename, poSrcDS, FALSE, 
                                    NULL, NULL, NULL );
    if( poDstDS != NULL )
        delete poDstDS;

In C:

    GDALDatasetH hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
    GDALDatasetH hDstDS;

    hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE, 
                             NULL, NULL, NULL );
    if( hDstDS != NULL )
        GDALClose( hDstDS );

In Python:

    src_ds = gdal.Open( src_filename )
    dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 )

CreateCopy()返回一个可写入的dataset，并且返回的dataset最终需要用户 自己关闭（和delete）以保证数据被真正地写入磁盘 （dataset本身可能有缓冲）。 参数FALSE表示当转换到输出格式时遇到不匹配或者丢失数据时，CreateCopy()宽大处理。 这主要是因为输 出格式可能不支持输入的数据类型，或者是不支持写操作。

一个更复杂的处理方式是指定某些选项，并且用预定义的回调函数获得进度。

In 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 );
    if( poDstDS != NULL )
        delete poDstDS;

In 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 );
    if( hDstDS != NULL )
        GDALClose( hDstDS );

In Python:

    src_ds = gdal.Open( src_filename )
    dst_ds = driver.CreateCopy( dst_filename, src_ds, 0, 
                                [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )

使用Create()
如果你不是简单地复制一个文件的话，就可能需要使用GDALDriver::Create() 来创建文件。Create()的参数列表和CreateCopy()相似，但是需要明确指定影象的大小、 波段数以及波段数据类型。
In C++:

    GDALDataset *poDstDS;       
    char **papszOptions = NULL;

    poDstDS = poDriver->Create( pszDstFilename, 512, 512, 1, GDT_Byte, 
                                papszOptions );

In C:

    GDALDatasetH hDstDS;        
    char **papszOptions = NULL;

    hDstDS = GDALCreate( hDriver, pszDstFilename, 512, 512, 1, GDT_Byte, 
                         papszOptions );

In Python:

    dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )

当dataset被正确地创建之后，特定的元数据和光栅数据都要被写到文件中。 这些操作一般需要依赖用户的具体选择，下边的代码是一个简单示例。

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

    delete poDstDS;

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

    GDALClose( hDstDS );

In Python:

    import Numeric, osr

    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 = Numeric.zeros( (512, 512) )    
    dst_ds.GetRasterBand(1).WriteArray( raster )

 

在打开GDAL所支持的光栅数据之前需要注册驱动。这里的驱动是针对GDAL支持的所有 数据格式。通常可以通过调用 GDALAllRegister() 函数来注册所有已知的驱动，同时 也包含那些用 GDALDriverManager::AutoLoadDrivers() 从.so文件中自动装载驱动。 如果程序需要对某些驱动做限制，可以参考 gdalallregister.cpp 代码。
当驱动被注册之后，我们就可以用 GDALOpen() 函数来打开一个数据集。打开的方式 可以是 GA_ReadOnly 或者 GA_Update。

In C++:

#include "gdal_priv.h"

int main()
{
    GDALDataset  *poDataset;

    GDALAllRegister();

    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
    if( poDataset == NULL )
    {
        ...;
    }

In C:

#include "gdal.h"

int main()
{
    GDALDatasetH  hDataset;

    GDALAllRegister();

    hDataset = GDALOpen( pszFilename, GA_ReadOnly );
    if( hDataset == NULL )
    {
        ...;
    }

In Python:

    import gdal
    from gdalconst import *

    dataset = gdal.Open( filename, GA_ReadOnly )
    if dataset is None:
        ...

如果 GDALOpen() 函数返回NULL则表示打开失败，同时 CPLError() 函数产生相应的错误信息。 如果您需要对错误进行处理可以参考 CPLError() 相关文档。通常情况下，所有的 GDAL函数都通过CPLError()报告错误。另外需要注意的是pszFilename并不一定对应一个 实际的文件名（当然也可以就是一个文件名）。它的具体解释由相应的驱动程序负责。 它可能是一个URL，或者是文件名以后后面带有许多用于控制打开方式的参数。通常建议， 不要在打开文件的选择对话框中对文件的类型做太多的限制。

获取Dataset信息
如果GDAL数据模型一节所描述的，一个GDALDataset包含了光栅数据的一系列的波段信息。 同时它还包含元数据、一个坐标系统、投影类型、光栅的大小以及其他许多信息。

    adfGeoTransform[0] /* 左上角 x */
    adfGeoTransform[1] /* 东西方向一个像素对应的距离 */
    adfGeoTransform[2] /* 旋转, 0表示上面为北方 */
    adfGeoTransform[3] /* 左上角 y */
    adfGeoTransform[4] /* 旋转, 0表示上面为北方 */
    adfGeoTransform[5] /* 南北方向一个像素对应的距离 */

如果需要输出dataset的基本信息，可以这样：

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

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

In Python:

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

获取一个光栅波段
现在，我们可以通过GDAL获取光栅的一个波段。同样每个波段含有元数据、块大小、 颜色表以前其他一些信息。下面的代码从dataset获取一个GDALRasterBand对象， 并且显示它的一些信息。
In 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() );

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

In Python:

        band = dataset.GetRasterBand(1)

        print 'Band Type=',gdal.GetDataTypeName(band.DataType)

        min = band.GetMinimum()
        max = band.GetMaximum()
        if min is not None and max is not None:
            (min,max) = ComputeRasterMinMax(1)
        print 'Min=%.3f, Max=%.3f' % (min,max)

        if band.GetOverviewCount() > 0:
            print 'Band has ', band.GetOverviewCount(), ' overviews.'

        if not band.GetRasterColorTable() is None:
            print 'Band has a color table with ', \
            band.GetRasterColorTable().GetCount(), ' entries.'

读光栅数据
GDAL有几种读光栅数据的方法，但是GDALRasterBand::RasterIO()是最常用的一种。 该函数可以自动转换数据类型、采样以及裁剪。下面的代码读光栅的第1行数据， 同时转换为float保存到缓冲。
In 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 );

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

In Python:

        scanline = band.ReadRaster( 0, 0, band.XSize, 1, \
                                     band.XSize, 1, GDT_Float32 )

返回的是一个string，包含了xsize*4大小的二进制数据，是float类型指针。 可以使用python的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 )

RasterIO()可以通过指定eRWFlag参数来确定是读/写数据（GF_Read或GF_Write）。 参数nXOff/nYOff/nXSize/nYSize描述了要读的影象范围（或者是写）。同时它也可以 自动处理边界等特殊情况。

参数pData指定读/写对应的缓冲。缓冲的类型必须是eBufType中定义的， 例如GDT_Float32、GDT_Byte等。RasterIO ()会自动转换缓冲和波段的类型， 使它们一致。当数据向下转换时，或者是数据超出转换后的数据类型可以表示的范围时， 将会用最接近的数据来代替。例如一个 16位的整数被转换为GDT_Byte时，所有大于255的 值都会用255代替（数据并不会被缩放）。

参数nBufXSize和nBufYSize描述了缓冲的大小。当时读写是是全部数据时， 该值和影象的大小相同。当需要对影象抽样的时候，缓冲也可以比真实的影象小。 因此，利用RasterIO()实现预览功能是很方便的。

参数nPixelSpace和nLineSpace通常被设置为0。当然，也可以使用他们来控制内存中的数据。 关闭Dataset

需要强调的一点是：GDALRasterBand对象属于相应的dataset，用户不能私自delete 任何GDALRasterBand对象。GDALDataset可以用GDALClose()关闭数据，或者是直接 delete GDALDataset对象。关闭GDALDataset的时候会进行相关的清除操作和刷新一些写操作。

创建文件的技巧
如果相应格式的驱动支持写操作的话，则可以创建文件。GDAL有两函数可以创建文件： CreateCopy()和Create()。 CreateCopy()函数直接从参数给定的数据集复制数据。 Create()函数则需要用户明确地写入各种数据（元数据、光栅数据等）。所有支持创建 的格式驱动都支持CreateCopy()函数，但是并不一定支持Create()函数。
为了确定数据格式是否支持Create或CreateCopy，可以检查驱动对象中的DCAP_CREATE 和DCAP_CREATECOPY元数据。在使用GetDriverByName()函数之前确保GDALAllRegister() 已经被调用过。

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

In C:

#include "cpl_string.h"
...
    const char *pszFormat = "GTiff";
    GDALDriver 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 );

In Python:

    format = "GTiff"
    driver = gdal.GetDriverByName( format )
    metadata = driver.GetMetadata()
    if metadata.has_key(gdal.DCAP_CREATE) \
       and metadata[gdal.DCAP_CREATE] == 'YES':
        print 'Driver %s supports Create() method.' % format
    if metadata.has_key(gdal.DCAP_CREATECOPY) \
       and metadata[gdal.DCAP_CREATECOPY] == 'YES':
        print 'Driver %s supports CreateCopy() method.' % format

我们可以看出有些格式不支持Create()或CreateCopy()调用。

使用CreateCopy()
GDALDriver::CreateCopy()函数使用比较简单，并且原先数据中的所有信息都被正确 的设置。函数还可以 指定某些可的选择参数，也通过一个回调函数来获得数据复制的 进展情况。下面的程序用默认的方式copy一个pszSrcFilename文件，保存 为 pszDstFilename 文件。
In C++:

    GDALDataset *poSrcDS = 
       (GDALDataset *) GDALOpen( pszSrcFilename, GA_ReadOnly );
    GDALDataset *poDstDS;

    poDstDS = poDriver->CreateCopy( pszDstFilename, poSrcDS, FALSE, 
                                    NULL, NULL, NULL );
    if( poDstDS != NULL )
        delete poDstDS;

In C:

    GDALDatasetH hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );
    GDALDatasetH hDstDS;

    hDstDS = GDALCreateCopy( hDriver, pszDstFilename, hSrcDS, FALSE, 
                             NULL, NULL, NULL );
    if( hDstDS != NULL )
        GDALClose( hDstDS );

In Python:

    src_ds = gdal.Open( src_filename )
    dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 )

CreateCopy()返回一个可写入的dataset，并且返回的dataset最终需要用户 自己关闭（和delete）以保证数据被真正地写入磁盘 （dataset本身可能有缓冲）。 参数FALSE表示当转换到输出格式时遇到不匹配或者丢失数据时，CreateCopy()宽大处理。 这主要是因为输 出格式可能不支持输入的数据类型，或者是不支持写操作。

一个更复杂的处理方式是指定某些选项，并且用预定义的回调函数获得进度。

In 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 );
    if( poDstDS != NULL )
        delete poDstDS;

In 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 );
    if( hDstDS != NULL )
        GDALClose( hDstDS );

In Python:

    src_ds = gdal.Open( src_filename )
    dst_ds = driver.CreateCopy( dst_filename, src_ds, 0, 
                                [ 'TILED=YES', 'COMPRESS=PACKBITS' ] )

使用Create()
如果你不是简单地复制一个文件的话，就可能需要使用GDALDriver::Create() 来创建文件。Create()的参数列表和CreateCopy()相似，但是需要明确指定影象的大小、 波段数以及波段数据类型。
In C++:

    GDALDataset *poDstDS;       
    char **papszOptions = NULL;

    poDstDS = poDriver->Create( pszDstFilename, 512, 512, 1, GDT_Byte, 
                                papszOptions );

In C:

    GDALDatasetH hDstDS;        
    char **papszOptions = NULL;

    hDstDS = GDALCreate( hDriver, pszDstFilename, 512, 512, 1, GDT_Byte, 
                         papszOptions );

In Python:

    dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )

当dataset被正确地创建之后，特定的元数据和光栅数据都要被写到文件中。 这些操作一般需要依赖用户的具体选择，下边的代码是一个简单示例。

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

    delete poDstDS;

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

    GDALClose( hDstDS );

In Python:

    import Numeric, osr

    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 = Numeric.zeros( (512, 512) )    
    dst_ds.GetRasterBand(1).WriteArray( raster )