gdal_retileで航空写真の分割

前回とは逆に、1枚物の航空写真(GeoTiff)をピクセルサイズ指定で分割するには「gdal_retile」を使うと良いんではないかと

gdal_retile -ps 5000 4000 -targetDir [出力先ディレクトリ名] target.tif

こんな感じで指定すると、target.tifファイルを横5000ピクセル、縦4000ピクセルの画像に分割し、[出力先ディレクトリ]の下に、「target_x_y.tif(x,yは自動で数字が振られる)」というファイル名で自動的に 分割されたファイルが生成されます。

  • targetDirオプションは必須っぽい

ここでもまた、tifファイルをLZW圧縮(高圧縮)したいのであれば

gdal_retile -ps 5000 4000 -co "COMPRESS=LZW" -co "PREDICTOR=2" -targetDir [出力先ディレクトリ名] target.tif

こんな感じで良さそう

gdal_mergeで航空写真の連結

分割されている航空写真(GeoTiff)などを連結する場合は「gdal_merge」コマンドを使うと良さそう

ただつなぎ合わせるだけであれば、

gdal_merge -o output.tif target1.tif target2.tif target3.tif target4.tif .....

のような感じで

  • o オプションの後に出力ファイル名、その後に連結したいファイルを半角スペース区切りで羅列していきます

連結時にLZW圧縮(高圧縮)したいのであれば

gdal_merge -o output.tif -co "COMPRESS=LZW" -co "PREDICTOR=2" target1.tif target2.tif target3.tif target4.tif .....

のような感じで指定すればOK

ogr2ogrでCSV形式の座標値リストの投影法変換

WGS84のUTM座標系(EPSG:32654)のX,Y,Zの情報が記録されたCSVファイルのうちX,Yをの平面直角座標12系(JGD2000)(EPSG:2454)に変換するよ

使用するのはogr2ogr フォーマットの変換や投影法の変換をするためのコマンドだけど、CSV形式の座標値を投影法だけ変換してCSVで出力するなんて使い方も出来るんだね

元のCSVファイルの中身はというと

[test.csv]-----------------------------------
xx,yy,z
609496.8867,4863616.468,172.3973083
609544.3473,4863616.582,166.1200714
609520.7308,4863616.97,166.2806549
609498.2774,4863617.479,171.5943604

                                                                                        • -

こんな感じ。

CSVファイルの変換をするためには、CSVと同じフォルダにこんな感じでvrtファイルを作成してあげて、

[test.vrt]-----------------------------------

  
    test.csv
    wkbPoint
    EPSG:32654
    
  

                                                                                        • -

コマンドラインCSVファイルのあるフォルダにカレントを移動して

ogr2ogr -f "CSV" -t_srs "EPSG:2454" -lco GEOMETRY=AS_XY output test.vrt

としてやると、outputフォルダが作成され、その中に「test.csv」が作成される。
中を確認してみると

[test.csv]-----------------------------------
X,Y,xx,yy,z
9142.8601224795475,-9156.8444376394891,609496.8867,4863616.468,172.3973083
9190.3242993202675,-9157.4487906071499,609544.3473,4863616.582,166.1200714
9166.7127514967997,-9156.7033192338477,609520.7308,4863616.97,166.2806549
9144.2661791730316,-9155.8544475998297,609498.2774,4863617.479,171.5943604

                                                                                        • -

こんな感じで、変換後のX,Y座標、その後ろに返還前のX,Y,Z座標というCSVファイルが作成されたよ。

gdal_translateして生成したGeoTiffが元のファイルサイズより大きくなってしまう

【元の画像】
 「base.tif」 :Size:122907KByte EPSG:32654 (GeoTiff・LZW圧縮)

photoshop等で航空写真.tifを編集・保存するとGeoTiff情報が欠落してしまうので、gdal_translateで座標系の情報を付与するため、コマンドラインより

gdal_translate -a_srs "EPSG:32654" base.tif work.tif

を実行すると、生成される「work.tif」のサイズが「430329KByte」
元のファイルサイズの3〜4倍もの大きさに。

上のコマンドだとLZW圧縮されていないためファイルサイズが大きくなってしまう
そこでオプションを追加

gdal_translate -a_srs "EPSG:32654" -co "COMPRESS=LZW" base.tif work.tif

よしよし、これで問題なかろうと思っていたら、生成されたwork.tifのサイズが「226680KByte」
非圧縮の場合よりは小さくなってるが、元画像と比較するとまだデカい…

そんな場合は、下記のように「-co "PREDICTOR=2"」オプションを追加すると良さそう

gdal_translate -a_srs "EPSG:32654" -co "COMPRESS=LZW" -co "PREDICTOR=2" base.tif work.tif

生成された「work.tif」のサイズは「122907KByte」

うむ、良い感じ。

LZW圧縮には2種類「低圧縮(PREDICTOR=1)」「高圧縮(PREDICTOR=2)」があるみたい
LZW圧縮(COMPRESS=LZW)時のデフォルトは「PREDICTOR=1」(低圧縮)となるようです

GDALで航空写真の投影法を変換する

ちょいと航空写真の投影法の変換をする用が出来たので手順をメモメモ
今回は、WGS84のUTM座標系(EPSG:32654)の航空写真を、JGD2000(世界測地系)の平面直角座標の12系(EPSG:2454)に変換するよ。


元の画像に座標系の情報が入っているか確認。

 gdalinfo <画像ファイル名>

表示される情報にEPSGうんちゃらとか入ってる?大丈夫?

もし投影法の情報が入ってないなら、gdal_translateコマンドで投影法の情報を付与する

 gdal_translate -a_srs "EPSG:32654" <入力ファイル名> <出力ファイル名>

後は投影法が付与された画像に対して、gdalwarpコマンドで投影法を変換する

 gdalwarp -t_srs "EPSG:2454" <入力ファイル名> <出力ファイル名>

以上

今回は元画像に投影法の情報を付与してからgdalwarpで変換したけど、元画像に投影法の情報が入って無くても

 gdalwarp -s_srs "EPSG:32654" -t_srs "EPSG:2454" <入力ファイル名> <出力ファイル名>

みたいに元画像の投影法をgdalwarpのオプションで指定して、直接変換することも可能なご様子

引越しだ。

某ブログシステムが使い辛いので引っ越してきた。
主に自分用の覚書を書いておく。

GDAL系のツールとかを使うこともままあるが、歳のせいか、コマンドやら引数やら覚える傍から脳みそから零れ落ちる日々なので、後から思い出せるように書いておきます。