GDALを使って地理座標を与えた地図画像を投影変換する

前回のエントリでは，Quantum GISを使って地図画像に地理座標を与えた。地理座標が定義できているので，地図投影法は自由に変えることが可能となる。今回は，希望の投影法に変換する手順を説明したい。

GDALは，ラスター（画像のようなメッシュ）形式のGISデータを処理するライブラリである。Linuxだけでなく，Windowsでも動作する。Windowsへのインストールは，まずOSGeo4Wで，OSGeo4W Installerをダウンロードし，実行させる。そして，GDALのみをチェックするとインストールできる。スタート > すべてのプログラム > OSGeo4W > OSGeo4W Shellを選択すると，コマンドプロンプトが出てきて，このプロンプトを通してGDALのライブラリが実行できる。コマンドプロンプトの使いかたについては，こちらのサイトが参考になりそうだ。

前回のエントリで行った，地理座標を与えた画像データがあるフォルダへ移動する。そして，前回の画像ファイルだが，Quantum GISによって座標はインプットされているのだが，投影法が画像ファイル自身には埋め込まれていないことが判明した。これは，gdalinfoコマンドで確認することができた。

gdalinfo <入力ファイル名>

投影法が定義されている場合には，投影法に関する記述が出てくるのだが，無いのである。なので，まずこれを定義することからはじめよう。定義するには，gdal_translateコマンドを使う。このコマンドは，画像の形式を変換することも可能な便利なコマンドである。

投影法の定義するためには，いろいろなパラメータが必要だが，コードで記述する方法が採用されている。それが，EPSGコードである。座標系と投影法の組み合わせで一つのEPSGコードが決まっている。なお，国内のEPSGコードについては，こちらのサイトで一覧が示されている。また，Google Earthと同じ，等緯度経度(Plate Caree)の場合のEPSGは4326，Google MapsはSpherical MercatorでEPSGが900913(または3785)である。

前回のエントリで扱った地図投影のEPSGコードはどうやって調べればよいかというと，前回のエントリの一つ目の画像が参考になる。WMSレイヤーを表示させるときに，基盤地図情報25000では，複数の投影法で表示させることが可能なことを記述したが，そこにEPSGのコードが出ている。それによると，25000分の1地形図の投影法はUTM図法のゾーン54の北側だったため，EPSGコードは3100である。それをgdal_translateでは以下のようにコマンドする。

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

出力ファイルはtifとすると，Geotiffとなるため，今回は全て拡張子はtifで統一する。出力ファイルをgdalinfoコマンドで表示させると，ちゃんと投影法が定義されていることが分かるはずである。これで，投影変換を行う準備が整った。投影変換のコマンドは以下のようになる。

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

-t_srs "EPSG:xxxx"のxxxxは変換先のEPSGコードである。なお，等緯度経度の場合はxxxxは4326となる。

実は，先ほどgdal_translateで元の画像のEPSGを定義したが，一気に変換することも可能である。その場合，EPSG:3100からEPSG:4326への変換となり，コマンドは以下のようになる。

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

"-s_srs"の後のEPSGが，入力ファイル名のEPSGコードである。s_srsとt_srsのsとtの違いは，sourceとtargetの違いである。それでは，出力されたファイルを比較してみよう。まずは，EPSG:3100(UTM/54N/JGD2000)である。

次がEPSG:4326である。

等緯度経度(Plate Caree)であるため，ゆがみが大きくなる。そのため，等緯度経度である下の画像の方が，横に伸びたような画像となっているのが分かる。

以上で，Quantum GISで地図画像に地理座標を与え，GDALで任意の投影変換を行う流れまでできたことになる。せっかくなので，次のエントリでは，等緯度経度に投影変換した地図画像をGoogle Earthでぴったりと重ねることができることを示したいと思う。