GDAL学習ノート-グリッドデータセットの読み取りと作成

グリッドファイルを開く
GDALがサポートするグリッドデータを開く前にドライバを登録する必要があります.GDALAllRegister()関数は、GDALDriverManager::AutoLoadDrivers()関数で自動的にロードされるすべての既知のドライバを登録しようとします.soファイル.

    GDALAllRegister();
    GDALDataset *poDataset;
    const char* pszFileName = "C:\\Users\\jinun\\Desktop\\SCRH\\dataset\\sub-TM-Spot-GS.img";
    poDataset = (GDALDataset *)GDALOpen( pszFileName, GA_ReadOnly);
    if( poDataset == NULL){
        printf( "Open the file failed!");
    }
    else{
        printf( "Open succeed!
" );
    }

GDALOpen()関数がNULLを返すと、オープンに失敗し、失敗情報はCPLError()関数で取得できます.
データセット情報の取得
adfGeoTransform[]を使用して、いくつかの情報を表すことができます.

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) */

次のコードは、データセットの情報を取得します.

    double adfGeoTransform[6];
    printf( "Driver: %s / %s
",
            poDataset->GetDriver()->GetDescription(),
            poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ));
    printf( "Size is %dx%dx%d
",
            poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
            poDataset->GetRasterCount());
    if ( poDataset->GetProjectionRef() != NULL)
    {
        printf( "Projection is '%s'
",
                poDataset->GetProjectionRef() );
    }
    if (poDataset->GetGeoTransform( adfGeoTransform ) == CE_None){
        printf( "Origin = (%.6f,%.6f)
",
                adfGeoTransform[0], adfGeoTransform[3] );
        printf( "Pixel Size = (%.6f,%.6f)
",
                adfGeoTransform[1], adfGeoTransform[5] );
    }

グリッド帯域の取得
グリッド帯域の取得はGDALで行うことができ、メタデータ、ブロックサイズ、カラーテーブル、その他の異なる情報も取得できます.

    GDALRasterBand *poBand;
    int nBlockXSize, nBlockYSize;
    int nGotMin, nGotMax;
    double adfMinMax[2];
    poBand = poDataset->GetRasterBand( 1 );
    poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
    printf( "Block = %dx%d Type = %s, ColorInterp = %s
",
            nBlockXSize, nBlockYSize,
            GDALGetDataTypeName( poBand->GetRasterDataType()),
            GDALGetColorInterpretationName( poBand->GetColorInterpretation()) );
    adfMinMax[0] = poBand->GetMinimum( &nGotMin );
    adfMinMax[1] = poBand->GetMaximum( &nGotMax );
    if( !( nGotMin && nGotMax )){
        GDALComputeRasterMinMax( (GDALRasterBandH)poBand, TRUE, adfMinMax );
    }
    printf( "Min = %.3fd, Max = %.3f
",adfMinMax[0], adfMinMax[1] );

    if( poBand->GetOverviewCount() > 0 ){
        printf( "Band has %d overviews.
", poBand->GetOverviewCount() );
    }
    if( poBand->GetColorTable() != NULL){
        printf( "Band has a color table with %d entries.
",
                poBand->GetColorTable()->GetColorEntryCount() );
    }

グリッドデータの読み込み
グリッドデータを読み取る方法はいくつかありますが、最も一般的な方法は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 );

GDALRasterBand::RasterIO()の使用については、このブログを参考にすることができます.
ファイルの作成
GDALでサポートされているファイルフォーマットでは、フォーマットドライバが作成をサポートしている場合は、新しいファイルを作成できます.新しいファイルを作成するには、CreateCopy()とCreate()の2つの方法があります.1つ目の方法はCreateCopy()関数を呼び出し、コピーが必要なソースデータセットパラメータを渡すことです.2つ目の方法はCreate()関数を呼び出すことです.新しいファイルを作成するすべてのドライバはCreateCopyメソッドをサポートし、一部のドライバのみがCreateメソッドをサポートします.駆動がCreateCopyメソッドをサポートするかCreateメソッドをサポートするかを決定するには、フォーマット駆動オブジェクトのDCAP_を確認する必要があります.CREATEとDCAP_CREATECOPYメタデータ.GetDriverByName関数を呼び出す前にGDALAllRegister()でドライバ登録する必要があります.次のコードは、ドライバがサポートする方法をテストします.

    GDALAllRegister();
    /*
     *   ，          ADRG、ECRGTOC、GEORASTER、GTiff、HDF4、HDF5、netCDF、
     * NITF、NTv2、OGDI、PDF、PostGISRaster、Rasterlite、RPFTOC、RS2、WCS WMS。
     */
    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.
", pszFormat );
    if ( CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATECOPY, FALSE) )
        printf( "Driver %s supports CreateCopy() method.
", pszFormat );

CreateCopyメソッドの使用
この方法では、コピーするファイルと保存するファイルを指定する必要があります.

 const char* pszSrcFileName = "C:\\Users\\jinun\\SCRH\\dataset\\sub-TM-Spot-GS.img";
    const char* pszDstFileName = "C:\\Users\\jinun\\SCRH\\dataset\\dst.img";;

    GDALDataset *poSrcDS = (GDALDataset *) GDALOpen( pszSrcFileName, GA_ReadOnly );
    GDALDataset *poDstDS;
    poDstDS = poDriver->CreateCopy( pszDstFileName, poSrcDS, FALSE,
                                    NULL, NULL, NULL );
    if( poDstDS != NULL){
        GDALClose( (GDALDatasetH) poDstDS );
    }
    GDALClose( (GDALDatasetH) poSrcDS );

より多くのパラメータを入力し、コピーの進捗状況を表示する必要がある場合は、次のように設計できます.

    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){
            GDALClose( (GDALDatasetH) poDstDS );
        }
        GDALClose( (GDALDatasetH) poSrcDS );

Createメソッドの使用
既存のファイルを指定した場所にコピーするだけでなく、GDALDriver::Create()メソッドを使用して新しいファイルを作成できます.create()メソッドはパラメータをcreatecopy()と同様に受け入れますが、画像のサイズ、帯域数、帯域タイプを提供する必要があります.

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

データセットの作成に成功すると、適切なメタデータとグリッドデータ情報はすべて新しいファイルに書き込まなければなりません.

double adfGeoTranform[6] = { 444720, 30, 0, 3751320, 0, -30 };
    OGRSpatialReference oSRS;//    #include "ogr_spatialref.h"   
    char *pszSRS_WKT = NULL;
    GDALRasterBand *poBand;
    GByte abyRaster[512 * 512];
    poDstDS->SetGeoTransform( adfGeoTranform );

    oSRS.SetUTM( 11, TRUE );
    oSRS.SetWellKnownGeogCS( "NAD27" );
    oSRS.exportToWkt( &pszSRS_WKT );
    CPLFree( pszSRS_WKT );

    poBand = poDstDS->GetRasterBand(1);
    poBand->RasterIO( GF_Write, 0, 0, 512, 512,
                      abyRaster, 512, 512, GDT_Byte, 0, 0);
    GDALClose( (GDALDatasetH) poDstDS );