まずは以下の４つの公式チュートリアルを参照してください。
これらのチュートリアルで大概のことは解決するはずです。さらには以下の２つのドライバプラグインも見ておくと中身がわかってよいでしょう。
国内ではこちらのサイトが詳しいです（が筆者にはレベルが高すぎ、、、）。
以下の小技集は、筆者が行き会ったものを適当に追加していきます。
GDALを利用するプログラムでの事前準備
GDALはイロイロなデータをファイルから読み込みますが、それらが格納されているディレクトリは環境変数で設定されている必要があります。
必要な環境変数は以下のとおりです。
- GDAL_DATA
- PROJ_LIB
- GEOTIFF_CSV
プログラムの中から設定する場合はCPLSetConfigOption()関数を使用します
（が、どうも反映されているようなきがしない、、、）。
1: #include <cpl_conv.h> 2: 3: CPLSetConfigOption( "GDAL_DATA", "C:/OSGeo4W/share/gdal" ); 4: CPLSetConfigOption( "PROJ_LIB", "C:/OSGeo4W/share/proj" ); 5: CPLSetConfigOption( "GEOTIFF_CSV", "C:/OSGeo4W/share/epsg_csv" );
GCPをセットするには
GCPをサポートするファイルフォーマットは少ないですが、GeoTIFFはサポートしていますので、
GCPを利用する場合はGeoTIFFを使用することになるでしょう。
以下のようにして設定します。
1: 2: // GCP数 3: int nGCP; 4: 5: // GCPの初期化 6: GDAL_GCP *pasGCPs = (GDAL_GCP *)CPLMalloc( nGCP * sizeof(GDAL_GCP) ); 7: GDALInitGCPs( nGCP, pasGCPs ); 8: 9: // 各GCPに値をセット 10: pasGCPs[0].dfGCPPixel = 画像座標 U; 11: pasGCPs[0].dfGCPLine = 画像座標 v; 12: pasGCPs[0].dfGCPX = GCP地上座標 X; 13: pasGCPs[0].dfGCPY = GCP地上座標 Y; ・ ・ ・ 14: 15: // データセットにGCPを設定 16: // poDSはGDALDataset型ポインタ 17: // pszProjectionRefは空間参照情報のWKT形式 18: poDS->SetGCPs( nGCPs, pasGCPs, pszProjectionRef );
GCP計算結果をシミュレートするには
GCPによる変換式だけがほしいときは、画像にGCPをセットする必要はありません。
例として、上記のようにGCPを詰め込んだ後、以下のようにしてGCPによる変換式を取得することができます。
double adfGeoTransform[6]; GDALGCPsToGeoTransform( nGCP, pasGCPs, adfGeoTransform, TRUE );
adfGeoTransform[]はdouble型6個の配列と決まっており、アフィン変換の係数が格納されます。
GDALGCPsToGeoTransform()第４引数はTRUEを指定すると多少の誤差を許容して最小二乗解を得ることができますが、
FALSEを指定した場合は誤差が0.25pixel以上あると、関数がFALSEを返して失敗するようになります。
GCPを用いてWARPを実行するには
画像にGCPをセットせずに、得られたGCPを用いて直接WARPを実行する場合は以下のようにします。
1: 2: // GCPTransformerを作成 3: // 第１引数のnGCPsはGCPの数 4: // 第２引数のpGCPはGCPデータ 5: // 第３引数は変換に使用する多項式の次数 6: // 第４引数はFALSEで順変換、TRUEで逆変換 7: void *hTransformArg = GDALCreateGCPTransformer( nGCPs, pGCP, 2, FALSE ); 8: 9: // 上記のハンドルを使って変換後の画像の地上座標、ピクセル数、ライン数を求める 10: double adfGeoTransform[6]; 11: CPLErr err = GDALSuggestedWarpOutput 12: ( hSrcDS, GDALGCPTransform, hTransformArg, adfGeoTransform, &nPixels, &nLines ); 13: 14: // ピクセル数、ライン数を求めたら出力画像を作成 15: GDALDataset *poDstDS = poDriver->Create 16: ( "out.tif", nPixels, nLines, GDALGetRasterCount( hSrcDS ), GDT_Byte, NULL ); 17: 18: ・ 19: ・ 20: ・ 21: 22: // ワープオプションに以下のように設定 23: psWarpOptions->pTransformerArg = hTransformArg; 24: psWarpOptions->pfnTransformer = GDALGCPTransform; 25: 26: ・ 27: ・ 28: ・ 29: 30: // ワープ実行後GCP変換のハンドルを解放する 31: GDALDestroyGCPTransformer( psWarpOptions->pTransformerArg ); 32: GDALDestroyWarpOptions( psWarpOptions );
ワープAPIについては公式チュートリアルを参照して下さい。
自動対応点探索を利用するには
GDALのクラスリストをぼやっと眺めているとGDALSimpleSURFなるクラスが有ります。
この名称を見てピンとくる向きも多いかと思いますが、これはSURF（Speed-Up Robust Features）を用いた自動マッチングで利用されるクラスです。
ただし、これらのクラスをGDALの外から利用することはできません。
自動マッチングはGDALComputeMatchingPoints()関数を用いて行います。
使い方は非常にシンプルで、以下のようにします。
1: char **papszOptions = NULL; 2: int nGCPs; 3: 4: // SURFで利用するDOG（Difference Of Gaussian）画像の開始解像度段階、終了解像度段階を指定する 5: // ピラミッド画像の同様のもので、SIFT/SURFではなぜかオクターブと呼ぶ 6: papszOptions = CSLSetNameValue( papszOptions, "OCTAVE_START", "1" ); 7: papszOptions = CSLSetNameValue( papszOptions, "OCTAVE_END", "5" ); 8: 9: // ２枚の画像のマッチングを行う 10: // GCPは１枚目の画像のピクセル／ラインと２枚目の画像の対応する点の地上座標となる 11: GDAL_GCP *pGCP = GDALComputeMatchingPoints( hSrcDS, hDstDS, papszOptions, &nGCPs );
あとは得られたGCPを使ってWARPを実行すると一応マッチしたものが得られるはずです。
が、誤マッチングを自動的に削除したりはしないので結果はそれなりです。
動作もやや安定していないようなので過度な期待は禁物です。
ある程度安定した結果を求めたい場合は、現状ではOpenCVで対応点探索→GCPデータに変換、としたほうが良さそうです。
SURFについての詳細は
コンピュータビジョン最先端ガイド2 [CVIMチュートリアルシリーズ]
あたりを参照して下さい。
WARP時の内挿方法を指定するには
Warpオプションを指定するところにあります。
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); psWarpOptions->eResampleAlg = GRA_Bilinear;
デフォルトはNearestなので注意しましょう。
カラーテーブルを取得・設定するには
以下の例は単バンド画像に対するカラーテーブルを取得する例です。
// poDatasetはGDALDataset型ポインタ GDALRasterBand *poBand = poDataset->GetRasterBand( 1 ); GDALColorTable *pColorTable = poBand->GetColorTable()->Clone()
GDALRasterBand::GetColorTable()関数が返すカラーテーブルのポインタは、
GDALRasterBandオブジェクトが所有権を持っており、読み取り専用です。
カラーテーブルを編集したいときは、上記のように複製を作成する必要があります。
設定は以下のようにします。
poBand->SetColorTable( pColorTable );
Nodataを取得／設定するには
いろいろな場面でNodataは出てきます。
Nodataは各バンドに対して個別に設定します。
// poBandはGDALRasterBand型ポインタ // 取得 // Nodataが設定されいればbSuccessにTRUEが返される int bSuccess; double dNodata = poBandSrc->GetNoDataValue( &bSuccess ); // 設定 poBand->SetNoDataValue( -9999.0 );
WARP時にもNodataを設定することが可能ですが、こちらは若干ややこしいです。
GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); psWarpOptions->padfSrcNoDataImag = (double *)CPLMalloc( psWarpOptions->nBandCount*sizeof(double) ); psWarpOptions->padfSrcNoDataImag[0] = 0.0; psWarpOptions->padfSrcNoDataReal = (double *)CPLMalloc( psWarpOptions->nBandCount*sizeof(double) ); psWarpOptions->padfSrcNoDataReal[0] = -9999.0; psWarpOptions->padfDstNoDataImag = (double *)CPLMalloc( psWarpOptions->nBandCount*sizeof(double) ); psWarpOptions->padfDstNoDataImag[0] = 0.0; psWarpOptions->padfDstNoDataReal = (double *)CPLMalloc( psWarpOptions->nBandCount*sizeof(double) ); psWarpOptions->padfDstNoDataReal[0] = -9999.0;
上記のように、ワープ前のNodataとワープ後のNodataをそれぞれ指定します。
なぜか複素数に対応しています。虚数部を省略することができません。
ラスタバンドの統計値を取得するには
バンドごとの簡単な統計値を取得することができます。
double dMin, dMax;
poBand->GetStatistics( FALSE, TRUE, &dMin, &dMax, NULL, NULL );
第１引数は、オーバービューなどが存在してかつその値でもいい場合はTRUEを指定します。
第２引数はTRUEを指定すると画像を一度走査します。FALSEが指定できるケースとは、ファイル中にこれらの値がすでに記述されている場合です。
まれにそのようなフォーマットがあるらしいです。
第３引数以降は取得したい値で、順番に最小値、最大値、平均値、標準偏差になっています。
不要な値はNULLを指定することができます。
自前のプログレス関数を設定するには
処理が長くなりそうないくつかの関数やワープオプションにはプログレスコールバック機構が用意されています。
これはどのように使うかというと、コールバック関数を指定しておくと定期的にその関数が処理途中に呼び出されるので、
画面などに「今何％」と表示することができるというものです。
コールバック関数は以下の形式です。
int (* GDALProgressFunc)(double dfComplete, const char *pszMessage, void *pProgressArg)
第１引数は進捗割合（0.0～1.0）、第２引数はGDALからのメッセージがそれぞれ渡されます。
第３引数は任意のオブジェクトで、この後出てきます。
GDALにはあらかじめ出来合いの関数として、GDALDummyProgress()、GDALTermProgress()、GDALScaledProgress()が用意されています。
以下の例は、Qtのプログレスバーに進捗を表示する例です。
1: 2: 3: class MainDialog : public QDialog, private Ui::DialogMain 4: { 5: Q_OBJECT 6: 7: public: 8: MainDialog(void); 9: ~MainDialog(void); 10: 11: // progressBarに値をセットする関数 12: inline void setProgress( int nComplete ) { progressBar->setValue( nComplete ); } 13: 14: // ・ 15: // ・ 16: // ・ 17: }; 18: 19: ///////////////////////////////////// 20: 21: // コールバック 22: static int CPL_STDCALL MyProgress( double dfComplete, const char *pszMessage, void *pProgressArg ) 23: { 24: ((MainDialog *)pProgressArg)->setProgress( (int)( dfComplete * 100 ) ); 25: return TRUE; 26: } 27: 28: 29: // MainDialog内からCreateCopy関数でコールバック 30: // thisポインタはMainDialogオブジェクト 31: // このポインタがコールバックの第３引数にそのまま渡される 32: 33: GDALDataset *poDstDS = poDriver->CreateCopy( strOutputFName.toLocal8Bit(), poVDS, 34: FALSE, papszOptions, MyProgress, this ); 35: 36: // MainDialog内からワープでコールバック 37: // ワープオプションで指定する 38: 39: GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); 40: 41: psWarpOptions->pfnProgress = MyProgress; 42: psWarpOptions->pProgressArg = this;
中央子午線を取得するには
投影座標系の中央子午線を取得するには以下のようにします。
// pszProjectionRefには空間参照のWKTが入っている OGRSpatialReference srs( pszProjectionRef ); double dPrimeMerid = srs.GetPrimeMeridian();
地理座標系を指定した場合は0.0が帰ってきたと思います。
仮想データセットを利用する
gdal_translateツールで使われているこの手法は非常に汎用性の高い方法です。
簡単に説明すると、実際にファイルを作成せずに仮想データセットを作成し、
必要な情報をあれこれ詰め込んだ後でCreateCopy()関数で一気に書き出してしまうという方法です。
CreateCopy()のほかにも、GCPを詰め込んでワープするといった用途にも使用できます。
したがって、gdal_translateでGCP詰め込み => gdalwarpで変換、の一連の操作を一気に行うことができます。
以下の例は仮想データセットにあれこれ詰め込んでCreateCopy()で書き出す例で、
gdal_translateツールとほぼ同じ手法です。
1: #include <vrt/vrtdataset.h> 2: 3: 4: // 仮想データセット作成 5: int nWidth = 500, nHeight = 500; 6: VRTDataset *poVDS = (VRTDataset *)VRTCreate( nWidth, nHeight ); 7: 8: // 位置情報セット 9: double adfGeoTransform[6] = { 100.0, 1.0, 0.0, -100.0, 0.0, -1.0 }; 10: poVDS->SetGeoTransform( adfGeoTransform ); 11: 12: // 地理参照をセット 13: char *pszProjectionRef = NULL; 14: 15: OGRSpatialReference oSRS; 16: oSRS.importFromEPSG( 2451 ); 17: oSRS.exportToWkt( &pszProjectionRef ); 18: 19: poVDS->SetProjection( pszProjectionRef ); 20: 21: // １バンド追加 22: poVDS->AddBand( GDT_UInt16, NULL ); 23: 24: // 追加したバンドを取得 25: VRTSourcedRasterBand *poVRTBand = (VRTSourcedRasterBand *)poVDS->GetRasterBand( 1 ); 26: 27: // 既存のデータセット（poDataset）からバンド情報をコピー 28: poVRTBand->AddSimpleSource( poDataset->GetRasterBand( 1 ) ); 29: 30: // 追加したバンドにNodataを設定 31: poVRTBand->SetNoDataValue( -9999.0 ); 32: 33: // CreateCopy()で書き出し 34: GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName( "GTiff" ); 35: char **papszOptions = NULL; 36: GDALDataset *poDstDS = poDriver->CreateCopy( "out.tif", poVDS, 37: FALSE, papszOptions, GDALTermProgress, NULL ); 38: 39: if ( poDstDS ) 40: GDALClose( (GDALDatasetH)poDstDS ); 41: 42: // 仮想データセットを破棄 43: GDALClose( (GDALDatasetH)poVDS ); 44:
“EPSG:4612″を解析するには
どういうことかというと、たとえばgdal_translateツールにある-a_srsオプションなんかは上記の形式を受け付けますが、
基本的にGDALが空間参照を取り扱うときはWKT形式です。したがってこのような空間参照の指定に対してWKTに変換することが必要になってきます。
これを行ってくれるのがOGRSpatialReference::SetFromUserInput()関数です。
見事にぴったりなものが用意されていますが、GDALは意外とこのような小技的関数が充実しています。
OGRSpatialReference srs; OGRErr err = srs.SetFromUserInput( "EPSG:4612" ); srs.exportToWkt( &pszProjectionRef );
上記のerrには、正しく解釈できた場合はOGRERR_NONEが、
失敗した場合はそれ以外が返されます。
この関数はそのうち非常に重大な意味を持つことになるでしょう。。。
それはまだ公開はためらわれるところですが、いずれその時が来たらまた追加したいと思います。
Geometryオブジェクトを作成するには
OGRプログラム内で使用するGeometryオブジェクトを作成する方法はいろいろあります。
ここではいくつかの例を挙げます。
WKTテキストから作成
1: 2: #include <ogrsf_frmts.h> 3: #include <ogr_spatialref.h> 4: #include <iostream> 5: 6: int main( int argc, char *argv[] ) 7: { 8: char szWkt[] = "POLYGON((135.0 35.0, 136.0 35.0, 136.0 36.0, 135.0 36.0, 135.0 35.0))"; 9: char *pszTemp = szWkt; 10: OGRGeometry *pGeom; 11: OGRSpatialReference oSRS; 12: 13: oSRS.importFromEPSG( 4326 ); 14: OGRErr err = OGRGeometryFactory::createFromWkt( &pszTemp, &oSRS, &pGeom ); 15: if ( err != OGRERR_NONE ) 16: { 17: std::cerr << "Failed to create." << std::endl; 18: exit( 1 ); 19: } 20: 21: // 使用後は消去 22: OGRGeometryFactory::destroyGeometry( pGeom ); 23: 24: return 0; 25: }
この方法は文字列操作がやりやすいPythonで使うと便利です。
srs.ImportFromEPSG(4612)
wkt = "POLYGON((" + west + " " + north + "," + west + " " + south + "," +
east + " " + south + "," + east + " " + north + "," + west + " " + north + "))"
rect = ogr.CreateGeometryFromWkt(wkt, srs)
OGRPointオブジェクトを作成
1: 2: #include <ogrsf_frmts.h> 3: #include <ogr_spatialref.h> 4: #include <iostream> 5: 6: int main( int argc, char *argv[] ) 7: { 8: OGRPoint *pGeom; 9: OGRSpatialReference oSRS; 10: oSRS.importFromEPSG( 4326 ); 11: 12: pGeom = dynamic_cast<OGRPoint *>( OGRGeometryFactory::createGeometry( wkbPoint ) ); 13: if ( !pGeom ) 14: { 15: std::cerr << "Failed to downcast" << std::endl; 16: exit( 1 ); 17: } 18: pGeom->setX( 135.0 ); 19: pGeom->setY( 35.0 ); 20: 21: OGRGeometryFactory::destroyGeometry( pGeom ); 22: 23: return 0; 24: }
OGRLineStringオブジェクトを作成
いくつか方法がありますが、点列が配列になっている場合は以下のようにすると効率がいいでしょう。
1: 2: #include <ogrsf_frmts.h> 3: #include <ogr_spatialref.h> 4: #include <iostream> 5: 6: int main( int argc, char *argv[] ) 7: { 8: OGRLineString *pGeom; 9: OGRSpatialReference oSRS; 10: oSRS.importFromEPSG( 4326 ); 11: 12: static double adfX[] = {135.0, 135.5, 136.0}; 13: static double adfY[] = {35.0, 35.5, 36.0}; 14: 15: pGeom = dynamic_cast<OGRLineString *>( OGRGeometryFactory::createGeometry( wkbLineString ) ); 16: if ( !pGeom ) 17: { 18: std::cerr << "Failed to downcast" << std::endl; 19: exit( 1 ); 20: } 21: 22: pGeom->setPoints( 3, adfX, adfY ); 23: 24: OGRGeometryFactory::destroyGeometry( pGeom ); 25: 26: return 0; 27: }
似たような方法でOGRPolygonオブジェクトも作成することができますが、ポリゴンの場合はもう少し面倒になるのでWKTから作ったほうが手っ取り早いような気がします。
ある図形と交差する要素を検索するには
えらく用途が限定された使い方ですが、例えば行政界ポリゴンを使ってある市町村が属する図郭を検索したりするときに使えます。
（筆者は図郭の検索ばかりやっているような気が、、、）。
図郭のファイルを開いている状態で、OGRLayerオブジェクトを取得したとして、以下のようにします。
1: OGRLayer *poLayer = poDS->GetLayerByName( "レイヤー名" ); 2: 3: // pGeomは対象図形 4: poLayer->SetSpatialFilter( pGeom ); 5: 6: // 検索された図形を走査 7: OGRFeature *poFeature = poLayer->GetNextFeature(); 8: while ( poFeature != NULL ) 9: { 10: OGRGeometryFactory::destroyGeometry( pGeom ); 11: poFeature = poLayer->GetNextFeature(); 12: }
この場合、検索される要素は対象図形と交差する要素です。
完全に含まれる要素とかを検索したい場合はOGRGeometry::Within()関数などを使う必要があります。
また、上記のコードを見るとわかるように、フィルタ適用後は検索対象要素以外は”なかったこと”扱いになります。
GDALのバージョンを取得するには
使用しているGDALのバージョンを取得する関数はGDALVersionInfo()で、引数に渡す文字列によって出力される文字列が異なります。
引数の種類は以下の4つです。
printf( "Version Number : %sn", GDALVersionInfo( "VERSION_NUM" ) );
printf( "Release Date : %sn", GDALVersionInfo( "RELEASE_DATE" ) );
printf( "Release Name : %sn", GDALVersionInfo( "RELEASE_NAME" ) );
printf( "--version : %sn", GDALVersionInfo( "arho" ) ); 上記以外の文字列の場合
実行結果は以下のようになります。
Version Number : 2020100 Release Date : 20170623 Release Name : 2.2.1 --version : GDAL 2.2.1, released 2017/06/23
spatialiteが4.1.2以降の場合
GDALから利用するspatialiteがバージョン4.1.2以降である場合、
プリプロセッサでSPATIALITE_412_OR_LATERを指定する必要があります。
これによっていくつかの追加関数が有効になります。
また、GDALをソースコードからコンパイルする場合でも指定しておく必要があります。