何かと便利なgdal/ogrユーティリティコマンドの小技集です。
とりあえず筆者が遭遇したものを適当に書いてあります。
詳細はgdalのドキュメントを参照してください。
ベクタデータの座標で、経緯度⇔座標値の変換をするには
ogr2ogrコマンドで、フォーマットの変換と同時に座標値も変換することができます。
以下の例は、平面直角座標9系のdgnファイルをGRS80の経緯度のkmlファイルに変換する例です。
> ogr2ogr -f KML -s_srs EPSG:2451 -t_srs EPSG:4019 out.kml in.dgn
Shapeファイルなら
Shapeファイルの場合、.prjファイルがあって正しく座標系の情報が記述されていれば、
-s_srsオプションを省略することができます。
gdalがサポートするファイル形式のサポート状況を調べるには
gdalがサポートするフォーマットの一覧は
gdalinfo --formatsと入力すると出てきます。
ここに表示される文字は他のコマンドでのフォーマット指定時の文字になります。
同様に、ogrinfo --formatsとするとogrがサポートするフォーマットの一覧が出力されます。
個別のフォーマットについてのサポート状況を出力するには
gdalinfo --format **と入力します（–formatsではない）。
**の部分にフォーマット名を指定します。
ogrにはこれに該当するものは無いようです。
このとき表示されるオプションリストは、そのフォーマットでファイルを書き出すときに指定可能なものです。
個別のコマンドでこれらは表示されないので、あらかじめ調べておく必要があります。
複数の隣接する画像を１つのファイルにまとめるには
例えば図郭ごとに画像があってそれらをひとつのファイルにまとめて全体を見たいという場合はgdal_merge.pyを使います。
当然、入力するファイルはGeoTIFF情報やTFWファイルなどで位置情報を持っている必要があります。
ワイルドカードは受け付けません。
互いの画像が重複する場合は、先に指定した画像ファイルが下になります。
コマンド例：
下の例は、４枚の画像を、XY範囲がそれぞれ-200~200、地上解像度1.0で合成する例です。
> gdal_merge -o out.tif -ps 1.0 1.0
-ul_lr -200.0 200.0 200.0 -200.0 in_1.tif in_2.tif in_3.tif in_4.tif
インデックスカラー対策（その１）
gdal_mergeでは、インデックスカラーの画像は普通にやるとグレースケールに勝手に変換されてしまいます。
回避策はpct2rgbコマンドで入力画像をRGBに変換後、gdal_mergeを実行し、
rgb2pctコマンドでインデックスカラーに戻すというのがありますが、
中間ファイルができてしまうのであまり上策とはいえません。
そこで、こんな場合はgdal_retileを使う手があります。
本来このコマンドは、複数の画像を指定したサイズに再分割するものですが、
サイズ指定で画像全体が収まるサイズを指定すれば１枚の画像が作られます。
あらかじめ、１枚にしたときの出力ファイルのサイズを調べて、
-psオプションでそのサイズを指定します。
また、このコマンドは再分割が目的なので、
出力指定はファイルではなく出力画像を格納するフォルダを指定します。
なので、あらかじめ出力用フォルダを作っておきます。
インデックスカラー対策（その２）
公式FAQの方法は以下のとおりです。
gdal_translate -of "VRT"で、どれか１枚の画像に対するVRTファイルを作成するgdal_mergeを実行（残念な画像ができる）- マージ後の画像に対して
gdal_translate -of "VRT"を実行し、VRTファイルを作成する - １．で作成したVRTファイルから<ColorTable>の部分をコピーし、３．のVRTファイルに貼り付ける
- ４．のVRTファイルを入力として
gdal_translateを実行する
ピラミッド画像を作成するには
ピラミッド画像を作るにはgdaladdoを使います。
特に凝った指定をしないならば、以下のようにコマンドを入力します。
> gdaladdo 入力ファイル名 2 4 8 16 ⏎
最後の数列はオリジナルサイズの何分の１のサイズにするかを指定します。
指定した数だけサブイメージが作られます。
出力フォーマットがGeoTIFFの場合は、ピラミッド画像は入力画像に埋め込まれます。
GeoTIFFの場合、埋め込まれるサブイメージは、元のイメージの圧縮形式は持ち越すようですが、
必ずタイル形式になり、タイルサイズも128×128で固定です。
これは、少なくともgdal-1.6ではgdaladdoではなくgdal内部でそういう指定にしてあるため、
気に入らない場合はgdal本体を書き換える必要がありそうです。
（該当箇所はfrmts/gtiff/geotiff.cppの2713〜2721行目）
いつからかわかりませんが、少なくともv1.9ではタイルサイズを変更することができるようになっています。
タイルサイズの変更はコマンドオプションではなく、環境変数GDAL_TIFF_OVR_BLOCKSIZEを設定して指定します。
指定しなければ128×128になります。
また、タイルサイズは２のべき乗数である必要があります。
例：
> set GDAL_TIFF_OVR_BLOCKSIZE=512 ⏎ > gdaladdo 入力ファイル名 2 4 8 16 ⏎
nodataを利用する
gdalwarpで画像を変換すると、万図にならずにデータがない箇所ができることがよくあります。
このようなピクセルにはnodataを指定しておかないと、
例えば上記のgdal_mergeでくっつけたりするときに邪魔になります。
nodataの指定は、入力側と出力側の２つがあり、
それぞれ-srcnodata、-dstnodataで指定します。
-srcnodataは、”nodataとして定義されていないピクセルをnodataとして扱う”
時に指定します。カラー画像など、チャンネルが複数ある画像の場合はWindows環境なら
-srcnodata "255 255 255"のようにダブルクォーテーションで囲って指定します。
-dstnodataは、出力画像で新たにnodataを指定するときに使用します。
その結果、nodataになるピクセルは、-srcnodataで指定したピクセルと、
変換の結果生じた余白になります。
-srcnodataを指定して、-dstnodataを指定しなかった場合、
nodataは持ち越されるわけではなく、出力画像ではピクセル値が0になります。
GeoTIFFの場合、nodataの定義はタグ番号42113で指定されています。
このタグはGDAL_NODATAというものらしいですが、これはGDALだけで利用可能なもので、
従ってtiffinfoで見ても番号だけの表示で、しかも複数チャンネルでの値を見ることができません。
だったらgdalinfoで見てみようということで見てみても、違う値が表示されたりします。
GDAL 1.8ではRGB各値が表示されるようになっています。
なお、IMGフォーマットの場合は、gdalinfoで正しく表示されるようです。
既にnodataが定義されている画像に対して、新たに-srcnodataでnodata
のピクセル値を追加し、-dstnodataでまた別のピクセル値を指定したりすると、なかなかすごいものができます。
こういうことをすると、出力フォーマットによって挙動が変わるようです。
v1.8から、NODATAの定義を除去することができるようになっています。
除去する場合は-a_nodata noneオプションを指定します。
マルチバンドの画像からフォルスカラー画像を作るには
こういうのはMultiSpecでやればいいのですが、gdal_translateでもできます。
-bオプションでチャンネルを選択します。
-a_nodataオプションでnodataを指定することもできます。
MultiSpecではピクセル値0のピクセルを背景にするかどうかの指定しかできないようです。
コマンド例は以下のような感じです。
> gdal_translate -ot "BYTE" -of "GTiff" -b 5 -b 2 -b 1
-a_srs "EPSG:2455" -a_nodata 0 "入力ファイル名" "出力ファイル名"
-bオプションはチャンネル数分指定します。
チャンネル番号は１スタートです。
gdal_translateのその他のオプションを使って、
かなり凝った指定もできます。
インデックスカラーのGeoTIFFをグレースケールに変換するには
見た目はグレースケールなのにフォーマット上はインデックスカラーになっていることがよくあります。
DIBフォーマットがこの方法でグレースケールを表現しています。
プログラムによってはインデックスカラーを受け付けないものもあるので、
このような場合はフォーマットを変換する必要がありますが、
レタッチソフトを通してしまうとGeoTIFF情報が失われるので困ります。
で、このような場合はgdal_translateで以下のように指定します。
> gdal_translate -of "GTiff" -co "PHOTOMETRIC=MINISBLACK"
"入力ファイル名" "出力ファイル名"
こうすればGeoTIFF情報は保存されます。
画像の地上解像度を下げるには
これもgdal_translateでできます。-outsizeオプションを使用します。
しかし、指定方法は地上解像度を指定するのではなく、変換後のピクセル／ライン数か、入力画像のピクセル数に対する比率で指定します。
したがって、例えば100x100pixで地上解像度が20cmの画像を地上解像度30cmにしたいという場合は以下のように指定する必要があります。
ピクセル／ライン数で指定の場合
> gdal_translate -outsize 67 67 "入力ファイル名" "出力ファイル名"
比率指定の場合
> gdal_translate -outsize 66.667% 66.667% "入力ファイル名" "出力ファイル名"
このように、割り切れない場合は端数が出てしまいます。
きっちり30cmにしたい、という場合は-ullrオプションで４隅範囲を指定して微調整するなどするしかない（？やや自信なし）と思われます。
上記のように縦横で別々の数値を指定するため、出力画像の縦横比は同じでなくてもかまいません
（わざわざそんなことをするケースもあまりないかとは思いますが）。
しかし、同時に縦横の地上解像度も別個に調整されるので、GISソフト上では画像の範囲は変わりません。
XY(Z)テキストデータを取り扱うには
測点データをとりあえずテキストファイルにした場合など、XY(Z)並びの点群データにはよく遭遇します。
このようなデータをShapeファイルにするには一度QGISを通せばできますが、
同時に座標変換をしたい場合などはできれば直接取り扱いたいところです。
ogrではCSV形式を扱うことができるので、これを利用します。
まず、テキストデータを何とかCSV形式にしてヘッダ行を追加し、以下のようにします。
x,y,z,Name 12.34,56.78,9.10,PT1 23.45,67.89,1.01,PT2 ・ ・ ・
自信はありませんが、”X”、”Y”の大文字は出力ファイルで使用される可能性があるので、
小文字にしておいたほうが無難です。
上記のファイルを「test.csv」というファイルに保存したとします。
さらに、VRTファイルというものを用意します。
gdal/ogrドキュメントのCSVフォーマットのページ
に説明がありますが、上記のファイルに対応したVRTファイルの内容は以下のようにします。
<OGRVRTDataSource>
<OGRVRTLayer name="test">
<SrcDataSource>test.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>WGS84</LayerSRS>
<GeometryField encoding="PointFromColumns" x="x" y="y"/>
</OGRVRTLayer>
</OGRVRTDataSource>
SrcDataSourceタグにファイル名を指定します。
また、GeometryFieldタグの”x=”、”y=”にそれぞれX座標とY座標が記述されている列名を指定します。
LayerSRSは、筆者は指定の仕方が良くわからないのでとりあえずこうしています。
わからなければコマンドの-s_srsオプションなどでSRSを指定すればいいので、まあ問題ないはずです。
例として、ogr2ogrコマンドでCSV点群ファイルをShapeファイルに変換しつつ、
平面直角座標６系からWGS84経緯度に変換するコマンドは以下のようにします。
> ogr2ogr -f "ESRI Shapefile" -s_srs EPSG:2448 -t_srs EPSG:4326 -lco GEOMETRY=AS_XY conv.shp test.vrt
入力ファイル名は.csvではなく.vrtを指定します。
CSVからCSVに変換は？
単に座標変換だけを行いたい場合はCSVからCSVに変換したいところです。
この場合は上記のコマンドで-f CSVとすればいいです。
ただし、gdal 1.6からのようです。
実行すると出力ファイル名で指定した名前のフォルダが作成され、その中にcsvファイルが作成されるので、
出力ファイル名ではなく出力フォルダ名をコマンドに入力します。
出力ファイルは、変換後の経緯度座標が先頭の２列にかきこまれ、その後に入力ファイルが並ぶ形になります。
フィールドの型を指定するには
CSVでは何も指定しないと文字列型として取り扱われてしまうようです。
整数型等の型として取り扱いたい場合は明示する必要があるようです。
この場合は<Field>タグで以下のように指定します。
<OGRVRTDataSource> <OGRVRTLayer name="points"> <SrcDataSource>points.csv</SrcDataSource> <SrcLayer>points</SrcLayer> <GeometryType>wkbPoint</GeometryType> <LayerSRS>WGS84</LayerSRS> <GeometryField encoding="PointFromColumns" x="x" y="y"/> <Field name="cat" type="Integer" /> <Field name="class" type="Integer" /> </OGRVRTLayer> </OGRVRTDataSource>
GDAL 2.1以降の場合
VRTを作成せずにコマンドラインオプションでXYZのフィールドを指定することができるようになりました。
-oo X_POSSIBLE_NAMES=オプションで以下のように指定します。
> ogr2ogr -f "ESRI Shapefile" -oo X_POSSIBLE_NAMES=X -oo Y_POSSIBLE_NAMES=Y out.shp input.csv
ヘッダ行がある場合は上記のようにフィールド名を指定します。
ヘッダ行がない場合は、１列目は’field_1’、２列目は’field_2’というように指定します。
既存のTIFF画像をGeoTIFFにするには
gdal_translateを利用する方法が２通りあります。
gdal_translate -a_ullrを利用する方法
このオプションを利用する場合は画像の左上座標と右下座標を指定します。
以下の例は平面直角座標２系のGeoTIFFを作成する例です。
> gdal_translate -of "GTiff" -a_srs EPSG:2444 -a_ullr 左上X座標 左上Y座標 右下X座標 右下Y座標 input.tif output.tif
当たり前ですが、右下座標を間違うと画像の縦横比率が変わってしまうので注意です。
tfwを作成して利用する方法
左上座標と地上解像度がわかっているような場合では、それらから右下の座標を計算して上記コマンドを入力してもいいですが、
tfwファイルを作成する方法も使えます。
tfwファイルの書式はこちらを参照してください。
tfwを作成したら、上記のコマンドから-a_ullrオプションを除いたコマンドを入力するとGeoTIFFが作成されます。
なお、上記の方法はいずれも画像が北上の場合のみ利用可能です。
回転も含めて変換したい場合はALOS L1BプロダクトをGeoTIFFに変換するにはを参照してください。
３つの画像からカラー画像を合成するには
GDALのツールを使ってカラー合成を行うには、VRTファイルを利用します。
VRTファイルはいろいろな使い方ができるXMLファイルですが、手書きでこれを作るのはまあ大変です。
ですが、カラー合成用のVRTファイルを作る場合はgdalbuildvrtコマンドでVRTファイルを作ることができます。
カラー合成用VRTファイルを作成する場合は、-separateオプションを使用します（v1.7から）。
以下のようにコマンドを入力します。
> gdalbuildvrt -separate "出力ファイル名"
"Red画像ファイル名" "Green画像ファイル名" "Blue画像ファイル名"
出来上がったVRTファイルは以下のような内容になっています。
例：
<VRTDataset rasterXSize="8401" rasterYSize="7261">
<SRS> SRS定義文 </SRS>
<GeoTransform> 地理座標変換係数 </GeoTransform>
<VRTRasterBand dataType="Byte" band="1">
<SimpleSource>
<SourceFilename relativeToVRT="1"> Red画像ファイル名 </SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="8401" RasterYSize="7261" DataType="Byte" BlockXSize="8401" BlockYSize="1" />
<SrcRect xOff="0" yOff="0" xSize="8401" ySize="7261" />
<DstRect xOff="0" yOff="0" xSize="8401" ySize="7261" />
</SimpleSource>
</VRTRasterBand>
<VRTRasterBand dataType="Byte" band="2">
<SimpleSource>
<SourceFilename relativeToVRT="1"> Green画像ファイル名 </SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="8401" RasterYSize="7261" DataType="Byte" BlockXSize="8401" BlockYSize="1" />
<SrcRect xOff="0" yOff="0" xSize="8401" ySize="7261" />
<DstRect xOff="0" yOff="0" xSize="8401" ySize="7261" />
</SimpleSource>
</VRTRasterBand>
<VRTRasterBand dataType="Byte" band="3">
<SimpleSource>
<SourceFilename relativeToVRT="1"> Blue画像ファイル名 </SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="8401" RasterYSize="7261" DataType="Byte" BlockXSize="8401" BlockYSize="1" />
<SrcRect xOff="0" yOff="0" xSize="8401" ySize="7261" />
<DstRect xOff="0" yOff="0" xSize="8401" ySize="7261" />
</SimpleSource>
</VRTRasterBand>
</VRTDataset>
このVRTファイルをQGISなどで読み込むとカラー合成された状態で表示されるので、
VRTファイル対応のソフトで閲覧するだけならこのままでも良いですが、
カラー画像ファイルを作りたい場合はこの後gdal_translateコマンドで変換すると良いでしょう。
> gdal_translate -of "GTiff" "VRTファイル名" "出力ファイル名"
ピクセル配列はセパレートではなく、BIP配列で作成されます。
複数のベクタデータファイルをマージするには
ogr2ogrコマンドの-appendオプションをつけて出力ファイルに既存のファイルを指定すると、
既存のファイルに変換もとのファイルが追加されます。
なので、複数のベクタデータに対して上記のコマンドをファイル数回実行すればマージすることができます。
例として、フォルダ内のShapeファイルを全て結合したい場合は、Dosコマンドならば
> for %i in (*.shp) do ogr2ogr -f "ESRI Shapefile" -append ..merge.shp %i
等とすればよいでしょう。出力ファイルが存在しない場合は-appendオプションは無視されます。
うまく利用すればフォーマットや座標系を変換しつつ複数のベクタデータファイルを一つにまとめることができるでしょう。
Shapeファイルのエンコードを変換するには
v1.9からShapeファイルのレイヤーオプションでエンコードを指定できるようになっています。
ogr2ogrコマンドで、
> ogr2ogr -f "ESRI Shapefile" -lco "ENCODING=UTF-8" 出力ファイル名 入力ファイル名
ただ、残念ながら現時点ではShift-JISへの変換はできないようです。
ベクタデータのフィールド情報を出力するには
今一つ使い方がわかりにくいogrinfoコマンドでベクタデータのフィールド情報のみを知りたい場合は以下のようにします。
> ogrinfo -ro -so -al -fields=yes "入力データ"
-fieldsオプションはデフォルトでYESなので、
-fields=yesはなくてもいいです。
このコマンドを実行すると、出力の最後のほうにフィールド情報が列挙されます。
マルチチャンネル画像のバンド順を入れ替えるには
例として、RGB画像をBGRに入れ替えて出力したい場合は以下のようにします。
> gdal_translate -b 3 -b 2 -b 1 input.tif output.tif
出力画像のバンド順は-bオプションの指定順になります。
したがって上記の場合は入力画像のバンド３が出力画像のバンド１、
入力画像のバンド１が出力画像のバンド３へとそれぞれ割り当てられます。