GeoTiff 座標系変換と結合

前提

国土地理院の地理院地図やGoogle Mapで提供されている航空写真を用いずに、独自に入手した衛星画像や航空写真を用いて広範囲の画像を表示させようとした時、座標系の変換と画像の結合が必要となります。

一般に入手される航空写真などで提供されている画像の座標系はJGD2011など、Web上では使い難い座標系を用いている場合が多く見受けられます。
その為、多くの場合には、座標系を変換する必要が生じるのです。
本資料では、この作業を行うための環境について説明し、衛星画像や航空写真などを用いて、座標系変換を行ってから画像を結合する以下の手順について説明します。

  1. 環境準備
  2. 画像の座標系を確認する。
  3. 座標系を変換する。

では順番に進めましょう。

環境準備

環境準備と言いましても、インストール方法については割愛します。
既に、本サイトでも紹介をしているOSGeo4Wをベースにお話しを進めます。
ただし、衛星画像や航空写真を結合する作業を行うには、正直申し上げてWindows環境はあまり適している環境とは言えません。
その理由としては以下が考えられるでしょう。

  • コマンドライン(シェル)が使い難い
  • 元々メモリの消費が激しく、大量のメモリを消費する作業では使い勝手が悪い。

数枚の画像を扱う程度であれば、GISツール(QGISなど)を用いて、ちまちま行うという方法もあるでしょう。パッチ処理などでまとめて出来るよ!という方もいらっしゃるとは思いますが。。。無視!

ということで少し作業がやり易いように1つツールを入れて作業を行っています。
また、一部の操作に関しては、Linux上で作業を行うかも知れません。

環境としては、以下の環境を準備しました。

Windows 10
OSGeo4W (GDALが入っていればOK!)
 参考情報:https://trac.osgeo.org/osgeo4w/wiki/OSGeo4W_jp
 過去の参考情報:https://tech.godpress.net/?p=368
Gow
 参考情報:https://ja.wikipedia.org/wiki/Gow

実際の作業ではGDALを使用します。
GDALのみをインストールした環境を準備されている方は、その環境を用いてもOKです。(って、そんな環境造っている人が、このサイトを参考にされるとは思いませんが。。。)

画像の座標系を確認する

OSGeo4Wをインストールしている場合、QGISで確認しちゃおう!という方は、GUI上でプロパティ確認で終わっちゃいますね。
でも今回は、コマンドライン操作が基本となりますので、コマンドで対処します。

Gowをインストールしていると、Linuxと同じような感覚でコマンドラインの操作が行える様になります。

なんでコマンドラインで作業を行うかと言いますと、入手した画像が全て同じ座標系であれば、あまり大したことではないのですが、時折、異なる座標系の画像を提供されて結合する場合があり、全ての画像を同じ座標系に統一しないと画像を結合した時にズレが生じて隙間が発生するのです。
間違って座標変換を行わないためには、事前に元画像の座標系を確認しておく必要があるためです。
では、実際のコマンド操作に関して説明します。

以下のコマンドは1つのファイルに対して座標系の確認を行うためのコマンドです。

複数のファイルに対し一括して確認を行う場合は、対象ファイルがあるフォルダーへ移動してから、以下のコマンドを実行します。

「gfind」コマンドは、Linuxの「find」コマンドをGowが実装したコマンドになります。

最初の引数「.」は対象となるフォルダーを指定しています。
対象フォルダーへ移動せずに、直接フォルダーを指定しても構いません。

-name *.tif :対象となるファイル名を指定します。ここではワイルドカードを用いた指定を行いました。全てのtifファイルが対象となります。

「|」で次のxargsコマンドへ出力結果を引き渡しています。
xargsコマンドでは-n 1とすることで、ファイル名を一つづつ処理します。
gdalinfo <ターゲットファイル>がファイル数分処理されます。

最後に、grepで必要な行だけを取得しています。

座標系を変換する

今回は、以下のフォルダー構成を想定して作業を進めます。
<作業フォルダー>—<EPSG102617>
         |
                                   +<EPSG4326>
コマンド操作は、<作業フォルダー>で行います。
<EPSG102617>は元の画像が保存されているフォルダーです。
先程の確認で座標系がEPSG:9001となっていましたが、内部ではEPSG:102617が正式な座標系となるためこの様にしています。

<EPSG4326>は座標変換後の画像ファイルを出力するためのフォルダーです。

既にお分かりのことと思いますが、本資料では、EPSG:102617(9001)→EPSG:4326の座標系変換を行います。

少し長くなりますが、以下のコマンドを実行してみましょう。

gfindコマンドで対象ファイルを特定し、そのファイル名だけを抜き出し、gdalwarpコマンドへ引き渡しています。

gdalwarpコマンドでは、-s_srsで元の座標系(102617)を-t_srsで指定される座標系(4326)へ変換を行うことを指定しています。
-srcnodataでは、元画像の中にある不要なデータのカラーコードを指定しています。同様に-dstnodataで出力時のノーデータを指定しています。

-srcnodataや-dstnodataを指定しなかった場合、画像の繋ぎ目に隙間が生じ、黒く筋が入ってしまうと思います。
座標変換に伴う歪により生じた隙間が黒色で塗潰され、それが隙間となってなって見えている状態になります。

これを除外するために、-srcnodataと-dstnodataを用いて不要なデータ部分を見えなくしているのです

余談・・・画像結合に関して

これだけの画像を準備したのですから、画像結合の話を少しだけ。。。
実は、画像結合してタイル画像を作成する予定だったのですが、マシンの不調とその必要性が無くなってしまったので、余談として記載します。

画像の結合にはgdal_mergeを使用します。

コマンドはこんな感じです。

-a_nodata “255 255 255″を追加した方が良さそうに思います。
環境変数GDAL_CHACHEMAXの値を出来るだけ大きくすると良さそうです。
ちなみに、gdalbuildvrtで仮想ファイルを作成して作業を行う方が効率良いとも思います。

こんな感じでしょうか?
タイル化などにご興味のある方は、別途ご連絡ください。