GdalTutorial GDAL使用简介

原文：http://www.gdal.org/gdal_tutorial.html

转载地址：http://blog.csdn.net/tiancai76/article/details/3063337

翻译：柴树杉

打开文件

在打开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 )

    {

        ...;

    }

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

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

    adfGeoTransform[0] 

    adfGeoTransform[1] 

    adfGeoTransform[2] 

    adfGeoTransform[3] 

    adfGeoTransform[4] 

    adfGeoTransform[5] 

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

In C++:

    double        adfGeoTransform[6];



    printf( "驱动：%s/%s/n",

            poDataset->GetDriver()->GetDescription(), 

            poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );



    printf( "数据大小：%dx%dx%d/n", 

            poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),

            poDataset->GetRasterCount() );



    if( poDataset->GetProjectionRef()  != NULL )

        printf( "投影类型：'%s'/n", poDataset->GetProjectionRef() );



    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )

    {

        printf( "起点：(%.6f,%.6f)/n",

                adfGeoTransform[0], adfGeoTransform[3] );



        printf( "像素大小：(%.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],                [3] );



        printf( "Pixel Size = (%.6f,%.6f)/n",

                adfGeoTransform[1], adfGeoTransform[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 ) ) );

读光栅数据

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

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

我们可以看出有些格式不支持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 );

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

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

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