ラスターAPIチュートリアル

ファイルを開きます

GDALでサポートされているラスターデータストアを開く前に,ドライバを登録する必要があります.各サポートされているフォーマットに対してドライバがあります.通常,これは,すべての既知のドライバを登録しようとする GDALAllRegister() 関数で行われます.これには. GDALDriverManager::AutoLoadDrivers() を使って.soファイルから自動的にロードされるドライバも含まれます.一部のアプリケーションでドライバのセットを制限する必要がある場合は.gdalallregister.cppからのコードを確認すると役立つかもしれません。Pythonは.gdalモジュールがインポートされると自動的にGDALAllRegister()を呼び出します.

ドライバが登録されたら,アプリケーションは,データセットを開くために独立した GDALOpen() 関数を呼び出すべきです.データセットの名前とアクセスしたいもの(GA_ReadOnly または GA_Update)を渡します.

C++の場合:

#include "gdal_priv.h"

#include <errno.h>

int main(int argc, const char* argv[])
{
    if (argc != 2) {
        return EINVAL;
    }
    const char* pszFilename = argv[1];

    GDALDatasetUniquePtr poDataset;
    GDALAllRegister();
    const GDALAccess eAccess = GA_ReadOnly;
    poDataset = GDALDatasetUniquePtr(GDALDataset::FromHandle(GDALOpen( pszFilename, eAccess )));
    if( !poDataset )
    {
        ...; // handle error
    }
    return 0;
}

Cの場合:

#include "gdal.h"
#include "cpl_conv.h" /* for CPLMalloc() */

#include <errno.h>

int main(int argc, const char* argv[])
{
    if (argc != 2) {
        return EINVAL;
    }
    const char* pszFilename = argv[1];

    GDALDatasetH  hDataset;
    GDALAllRegister();
    const GDALAccess eAccess = GA_ReadOnly;
    hDataset = GDALOpen( pszFilename, eAccess );
    if( hDataset == NULL )
    {
        ...; // handle error
    }
    return 0;
}

Pythonの場合:

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

Note that if GDALOpen() returns NULL it means the open failed, and that an error messages will already have been emitted via CPLError(). If you want to control how errors are reported to the user review the CPLError() documentation. Generally speaking all of GDAL uses CPLError() for error reporting. Also, note that pszFilename need not actually be the name of a physical file (though it usually is). It's interpretation is driver dependent, and it might be an URL, a filename with additional parameters added at the end controlling the open or almost anything. Please try not to limit GDAL file selection dialogs to only selecting physical files.

データセット情報の取得

As described in the ラスターデータモデル, a GDALDataset contains a list of raster bands, all pertaining to the same area, and having the same resolution. It also has metadata, a coordinate system, a georeferencing transform, size of raster and various other information.

特に,一般的な"北向き"画像で回転やシアリングがない場合,ジオリファレンス変換 Geotransform Tutorial は以下の形式を取ります:

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を介したラスターデータへのアクセスは,1つのバンドずつ行われます.また,メタデータ,ブロックサイズ,カラーテーブルおよびその他の情報がバンドごとに利用可能です.以下のコードは,データセットから 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)

返されたスキャンラインは文字列型であり,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の設定に基づいて読み込みまたは書き込みに使用する同じRasterIO()呼び出しが使用されることに注意してください(いずれかGF_ReadまたはGF_Write).nXOff, nYOff, nXSize, nYSize引数は,ディスク上のラスターデータのウィンドウを読み取る(または書き込む)ために使用されます.タイルの境界にある必要はありませんが,アクセスが効率的であるかもしれません.

pDataは,データが読み込まれるまたは書き込まれるメモリバッファです.実際の型は,eBufTypeとして渡されるものである必要があります.例えば,GDT_Float32またはGDT_Byteなどです. RasterIO() 呼び出しは,バッファのデータ型とバンドのデータ型の間の変換を処理します.浮動小数点データを整数に変換する場合,RasterIO()は切り捨てを行い,出力の許容範囲外のソース値を変換する場合,最も近い許容値が使用されます.これは,例えば,GDT_Byteバッファに読み込まれた16ビットデータは,255を超えるすべての値を255にマップし,データはスケーリングされません!

nBufXSizeとnBufYSizeの値は,バッファのサイズを示します.データを完全な解像度で読み込む場合,これはウィンドウサイズと同じになります.ただし,低解像度の概要を読み込む場合,これはディスク上のウィンドウよりも小さく設定することができます.この場合,RasterIO()は,適切な場合には概要を利用して,IOを効率的に行います.

通常はゼロであり,デフォルト値を使用することを示すnPixelSpaceとnLineSpaceです.ただし,他のピクセル間データを含むバッファに読み込むことを許可するために,メモリデータバッファへのアクセスを制御するために使用することができます.

データセットを閉じる

GDALRasterBand オブジェクトはデータセットによって所有されており,決してC++のdelete演算子で破棄してはいけません. GDALDataset は, GDALClose() を呼び出すことで閉じることができます(Windowsユーザーは,モジュール境界を越えてメモリを割り当ておよび解放する際の既知の問題のため,GDALDatasetに対してdelete演算子を使用することは推奨されません.FAQの関連トピックを参照してください).GDALCloseを呼び出すと,適切なクリーンアップと,保留中の書き込みのフラッシュが行われます.人気のあるフォーマットで開いたデータセットを更新モードで開いた場合にGDALCloseを呼び出し忘れると,その後開くことができなくなる可能性が高いです.

ファイルを作成するためのテクニック

GDALでサポートされているフォーマットで新しいファイルを作成することができます.フォーマットドライバが作成をサポートしている場合,ファイルを作成するための2つの一般的なテクニックがあります. CreateCopy() と Create() です. CreateCopyメソッドは,フォーマットドライバ上でCreateCopy()メソッドを呼び出し,コピーする必要があるソースデータセットを渡すことを含みます. Createメソッドは,ドライバ上でCreate()メソッドを呼び出し,すべてのメタデータとラスターデータを明示的に別々の呼び出しで書き込むことを含みます.新しいファイルを作成することをサポートするすべてのドライバは,CreateCopy()メソッドをサポートしていますが,Create()メソッドをサポートしているのはわずかです.

特定のフォーマットがCreateまたはCreateCopyをサポートしているかどうかを判断するには,フォーマットドライバオブジェクトのDCAP_CREATEおよびDCAP_CREATECOPYメタデータを確認することが可能です. GDALDriverManager::GetDriverByName() を呼び出す前に, GDALAllRegister() が呼び出されていることを確認してください.この例では,ドライバを取得し,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()の使用

The GDALDriver::CreateCopy() method can be used fairly simply as most information is collected from the source dataset. However, it includes options for passing format specific creation options, and for reporting progress to the user as a long dataset copy takes place. A simple copy from the a file named pszSrcFilename, to a new file named pszDstFilename using default options on a format whose driver was previously fetched might look like this:

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