GDAL Warp API チュートリアル (再投影、...)

概要

GDAL Warp API ( gdalwarper.h で宣言されています) は,アプリケーションが提供する地理的変換関数 (GDALTransformerFunc),さまざまなリサンプリングカーネル,およびさまざまなマスキングオプションを使用して,高性能の画像ワーピングサービスを提供します.メモリに保持できるよりもはるかに大きなファイルをワープできます.

このチュートリアルでは,ワープ API を使用したアプリケーションの実装方法を示しています. Warp API に対する C および Python のバインディングが不完全であることを前提としています. また ラスターデータモデル と一般的な GDAL API に精通していることを前提としています.

通常,アプリケーションは,使用するオプションを持つ GDALWarpOptions 構造体を初期化し,これらのオプションに基づいて GDALWarpOperation をインスタンス化し,その後,内部で GDALWarpKernel クラスを使用してワープオプションを実行するために GDALWarpOperation::ChunkAndWarpImage() メソッドを呼び出します.

単純な再投影ケース

最初に,適切な出力ファイルが既に存在し,エラーチェックが最小限であると仮定して,画像を再投影するための比較的単純な例を構築します.

#include "gdalwarper.h"

int main()
{
    GDALDatasetH  hSrcDS, hDstDS;

    // Open input and output files.
    GDALAllRegister();
    hSrcDS = GDALOpen( "in.tif", GA_ReadOnly );
    hDstDS = GDALOpen( "out.tif", GA_Update );

    // Setup warp options.
    GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();
    psWarpOptions->hSrcDS = hSrcDS;
    psWarpOptions->hDstDS = hDstDS;
    psWarpOptions->nBandCount = 1;
    psWarpOptions->panSrcBands =
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panSrcBands[0] = 1;
    psWarpOptions->panDstBands =
        (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
    psWarpOptions->panDstBands[0] = 1;
    psWarpOptions->pfnProgress = GDALTermProgress;

    // Establish reprojection transformer.
    psWarpOptions->pTransformerArg =
        GDALCreateGenImgProjTransformer( hSrcDS,
                                        GDALGetProjectionRef(hSrcDS),
                                        hDstDS,
                                        GDALGetProjectionRef(hDstDS),
                                        FALSE, 0.0, 1 );
    psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

    // Initialize and execute the warp operation.
    GDALWarpOperation oOperation;
    oOperation.Initialize( psWarpOptions );
    oOperation.ChunkAndWarpImage( 0, 0,
                                GDALGetRasterXSize( hDstDS ),
                                GDALGetRasterYSize( hDstDS ) );
    GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );
    GDALDestroyWarpOptions( psWarpOptions );
    GDALClose( hDstDS );
    GDALClose( hSrcDS );
    return 0;
}

この例では,既存の入力ファイルと出力ファイル (in.tif および out.tif) を開きます. GDALWarpOptions 構造体が割り当てられます (GDALCreateWarpOptions() は,たくさんのものに対して合理的なデフォルトを設定します,常にデフォルトのものに使用してください),そして,入力ファイルハンドルと出力ファイルハンドル,およびバンドリストが設定されます. panSrcBands と panDstBands リストはここで動的に割り当てられ, GDALDestroyWarpOptions() によって自動的に解放されます. シンプルなターミナル出力進行状況モニター (GDALTermProgress) がインストールされ,ユーザーに完了進行状況を報告します.

GDALCreateGenImgProjTransformer() is used to initialize the reprojection transformation between the source and destination images. We assume that they already have reasonable bounds and coordinate systems set. Use of GCPs is disabled.

オプション構造が準備できたら,それらを使用して GDALWarpOperation をインスタンス化し,実際に GDALWarpOperation::ChunkAndWarpImage() でワープを実行します. その後,変換器,ワープオプション,およびデータセットをクリーンアップします.

通常,ファイルを開いた後,再投影変換器を設定した後 (失敗した場合は NULL を返す),およびワープを初期化した後にエラーチェックが必要です.

その他のワーピングオプション

GDALWarpOptions 構造体には,ワーピング動作を制御するために設定できる項目がいくつか含まれています. 特に興味深いものはいくつかあります:

  • GDALWarpOptions::dfWarpMemoryLimit - Set the maximum amount of memory to be used by the GDALWarpOperation when selecting a size of image chunk to operate on. The value is in bytes, and the default is likely to be conservative (small). Increasing the chunk size can help substantially in some situations but care should be taken to ensure that this size, plus the GDAL cache size plus the working set of GDAL, your application and the operating system are less than the size of RAM or else excessive swapping is likely to interfere with performance. On a system with 256MB of RAM, a value of at least 64MB (roughly 64000000 bytes) is reasonable. Note that this value does not include the memory used by GDAL for low level block caching.

  • GDALWarpOptions::eResampleAlg - One of GRA_NearestNeighbour (the default, and fastest), GRA_Bilinear (2x2 bilinear resampling) or GRA_Cubic. The GRA_NearestNeighbour type should generally be used for thematic or color mapped images. The other resampling types may give better results for thematic images, especially when substantially changing resolution.

  • GDALWarpOptions::padfSrcNoDataReal - This array (one entry per band being processed) may be setup with a "nodata" value for each band if you wish to avoid having pixels of some background value copied to the destination image.

  • GDALWarpOptions::papszWarpOptions - This is a string list of NAME=VALUE options passed to the warper. See the GDALWarpOptions::papszWarpOptions docs for all options. Supported values include:

    • INIT_DEST=[value] または INIT_DEST=NO_DATA: このオプションは,宛先画像を指定された値 (すべてのバンドに対して) で初期化するか, padfDstNoDataReal/padfDstNoDataImag の NO_DATA 値で初期化することを示します. この値が設定されていない場合,宛先画像は読み込まれ,ソースワープが重ねられます.

    • WRITE_FLUSH=YES/NO: このオプションは,各チャンクが処理された後にデータをディスクにフラッシュすることを強制します. いくつかの場合,これにより,出力データの直列書き込みが確実になります. そうでない場合,入力バッファのためにデータのブロックが読み込まれるたびにデータのブロックがディスクに書き込まれ,ディスクの周りを多くの余分なシークし,IO スループットが低下します. この時点でのデフォルトは NO です.

出力ファイルの作成

前のケースでは,適切な出力ファイルが既に存在すると仮定されていました. ここでは,新しい座標系の適切な境界を持つ新しいファイルが作成されるケースを説明します. この操作は,ワープ API に特に関連していません. これは単に変換 API を使用しています.

#include "gdalwarper.h"
#include "ogr_spatialref.h"
...
GDALDriverH hDriver;
GDALDataType eDT;
GDALDatasetH hDstDS;
GDALDatasetH hSrcDS;

// Open the source file.
hSrcDS = GDALOpen( "in.tif", GA_ReadOnly );
CPLAssert( hSrcDS != NULL );

// Create output with same datatype as first input band.
eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS,1));

// Get output driver (GeoTIFF format)
hDriver = GDALGetDriverByName( "GTiff" );
CPLAssert( hDriver != NULL );

// Get Source coordinate system.
const char *pszSrcWKT, *pszDstWKT = NULL;
pszSrcWKT = GDALGetProjectionRef( hSrcDS );
CPLAssert( pszSrcWKT != NULL && strlen(pszSrcWKT) > 0 );

// Setup output coordinate system that is UTM 11 WGS84.
OGRSpatialReference oSRS;
oSRS.SetUTM( 11, TRUE );
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.exportToWkt( &pszDstWKT );

// Create a transformer that maps from source pixel/line coordinates
// to destination georeferenced coordinates (not destination
// pixel line).  We do that by omitting the destination dataset
// handle (setting it to NULL).
void *hTransformArg;
hTransformArg =
    GDALCreateGenImgProjTransformer( hSrcDS, pszSrcWKT, NULL, pszDstWKT,
                                     FALSE, 0, 1 );
CPLAssert( hTransformArg != NULL );

// Get approximate output georeferenced bounds and resolution for file.
double adfDstGeoTransform[6];
int nPixels=0, nLines=0;
CPLErr eErr;
eErr = GDALSuggestedWarpOutput( hSrcDS,
                                GDALGenImgProjTransform, hTransformArg,
                                adfDstGeoTransform, &nPixels, &nLines );
CPLAssert( eErr == CE_None );
GDALDestroyGenImgProjTransformer( hTransformArg );

// Create the output file.
hDstDS = GDALCreate( hDriver, "out.tif", nPixels, nLines,
                     GDALGetRasterCount(hSrcDS), eDT, NULL );
CPLAssert( hDstDS != NULL );

// Write out the projection definition.
GDALSetProjection( hDstDS, pszDstWKT );
GDALSetGeoTransform( hDstDS, adfDstGeoTransform );

// Copy the color table, if required.
GDALColorTableH hCT;
hCT = GDALGetRasterColorTable( GDALGetRasterBand(hSrcDS,1) );
if( hCT != NULL )
    GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
... proceed with warp as before ...

このロジックに関するいくつかの注意:

  • 出力座標がピクセル行座標ではなく,ジオリファレンスされた出力であるように出力座標を作成する必要があります. なぜなら,ソース画像の周りのピクセルを宛先ジオリファレンス座標にマッピングするために変換器を使用するからです.

  • The GDALSuggestedWarpOutput() function will return an adfDstGeoTransform, nPixels and nLines that describes an output image size and georeferenced extents that should hold all pixels from the source image. The resolution is intended to be comparable to the source, but the output pixels are always square regardless of the shape of input pixels.

  • ワーパーは,"ランダム"に書き込むことができる形式の出力ファイルを必要とします. これは,通常,Create() メソッドの実装を持つ非圧縮形式に制限されます (CreateCopy() とは対照的). 圧縮形式にワープするか,CreateCopy() スタイルの形式にワープするには,より適切な形式で画像の完全な一時コピーを作成し,それを望ましい最終形式に CreateCopy() する必要があります.

  • Warp API は,ピクセルのみをコピーします. すべてのカラーマップ,ジオリファレンシング,およびその他のメタデータは,アプリケーションによって宛先にコピーする必要があります.

パフォーマンスの最適化

ワープ API のパフォーマンスを最適化するために行うことができるいくつかのことがあります:

  • ワープ API チャンキングに使用可能なメモリ量を増やして,一度により大きなチャンクで操作できるようにします. これは GDALWarpOptions::dfWarpMemoryLimit パラメータです. 理論的には,操作されるチャンクサイズが大きいほど,より効率的な I/O 戦略と,より効率的な近似変換になります. ただし,ワープメモリと GDAL キャッシュの合計は RAM サイズよりも小さくなければならず,おそらく RAM サイズの 2/3 ほどでしょう.

  • GDAL キャッシュのメモリ量を増やします. これは,非常に大きな入力および出力画像をスキャンライン指向で処理する場合に特に重要です. 入力または出力のすべてのスキャンラインが各チャンクで交差するたびに再読み込みされる場合,パフォーマンスが大幅に低下する可能性があります. GDAL 内で使用可能なキャッシュのメモリ量を制御するには, GDALSetCacheMax() を使用します.

  • 変換される各ピクセルの正確な再投影の代わりに近似変換を使用します. このコードは,再投影変換に基づいて近似変換が作成される方法を示していますが,指定された誤差しきい値 (出力ピクセルの dfErrorThreshold) で行います.

hTransformArg =
    GDALCreateApproxTransformer( GDALGenImgProjTransform,
                                 hGenImgProjArg, dfErrorThreshold );
pfnTransformer = GDALApproxTransform;

  • 空の出力ファイルに書き込む場合,出力チャンクを固定値で初期化するために, GDALWarpOptions::papszWarpOptions の INIT_DEST オプションを使用します. これにより,不要な IO 処理が大幅に削減される可能性があります.

  • タイル形式の入力および出力形式を使用します. タイル形式を使用すると,追加の画像データに触れることなく,指定されたソースおよび宛先画像のチャンクにアクセスできます. 大きなスキャンライン指向のファイルは,多くの余分な IO が発生する可能性があります.

  • すべてのバンドを 1 回の呼び出しで処理します. これにより,変換計算が各バンドごとに実行される必要がなくなります.

  • Use the GDALWarpOperation::ChunkAndWarpMulti() method instead of GDALWarpOperation::ChunkAndWarpImage(). It uses a separate thread for the IO and the actual image warp operation allowing more effective use of CPU and IO bandwidth. For this to work GDAL needs to have been built with multi-threading support (default on Win32, default on Unix, for previous versions -with-threads was required in configure).

  • リサンプリングカーネルは,最も単純な最近傍補間から始まり,バイリニアおよび立方体でより複雑になります. 必要以上に複雑なリサンプリングカーネルを使用しないでください.

  • 一般的な特殊ケースに対して特別な簡略化されたロジックケースが使用されるように,エソテリックなマスキングオプションの使用を避けます. たとえば,8bit データでマスキングなしの最近傍補間は,一般的なケースと比較して高度に最適化されています.

その他のマスキングオプション

GDALWarpOptions には,入力および出力の有効性マスクおよび密度マスクのためのエソテリックなマスキング機能がいくつか含まれています. これらのうちいくつかはまだ実装されておらず,他のいくつかは実装されていますが,不十分にテストされています.バンドごとの有効性マスク以外の機能は,現時点ではこれらの機能を慎重に使用することが推奨されます.