GDAL API入門

pamxy發表於2013-04-14
API

轉自：http://www.osgeo.org.cn/l18n/gdal/gdal_tutorial.html

翻譯：柴樹杉(chaishushan@gmail.com) 
原文：http://www.gdal.org/gdal_tutorial.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 )

相關文章