基盤地図情報10mメッシュ(標高)をMapinfoでつかう

遅ればせながら基盤地図情報の10mメッシュ(標高)を使ってみました。
このデータ、せっかく無償で公開しているのに、ちょっと使いにくいです。地図センターなんて、基盤地図情報（数値標高モデル）変換サービスとかやっちゃっています。1メッシュあたり280円也。GTPOPOもSRTMも無償で公開されていたけれど、有償で変換する商売なんてあったかなぁ。
というわけで、なんとか無料でMapinfoで使いたい。その場合の課題としては。

基盤地図情報のメッシュ(標高)は緯度経度になっている
平面直角やUTMで使いたいがMapinfoはラスタの投影変換できない・・
いったんポイントにして投影変換した後、内挿してラスタを作り直すとかめんどうすぎる

実は会社でやればERDAS IMAGINEのインポータとかを使ってすぐできるのですが、それではつまらないし、IMAGINE持ってないとどうしようもなくなってしまうので、オープンソースのツールだけでなんとかMapinfoまで持って行く方向でやり方を調べました。

１）基盤地図情報DEMコンバータでxmlをbilに変換
ダウンロード・解凍したデータを、ソフト公開のページ (三匹のウリボウ)で公開されている基盤地図情報DEMコンバータを使って、bilに変換。このときに適当な範囲をマージしておくこともできます。ここではマージして、union.bilをいうファイル名で保存しました。赤木さま、便利なツールをありがとうございます。

２）gdalwarpで緯度経度から平面直角に投影変換
（１）でつくったbilファイルは緯度経度のままなので、これをMapinfoで表示したとしても南北につぶれてしまいます（実際にはMapinfo＋vertical mapperではbilのインポートはできませんが・・）。なので、ここでは緯度経度（世界測地系）→平面直角○系（世界測地系）のラスタ投影変換を行います。使うツールはgdalwarp。Windowsな人であればFWToolsを導入すればFWTools Shellから使うことができます。
↓みたいな感じ。

gdalwarp -s_srs EPSG:4612 -t_srs EPSG:2448 -tr 10 10 -r bilinear -dstnodata -9999 union.bil union_xy.tif

　-s_srs EPSG:4612はインプットファイルのEPSGコード（ここではunion.bil、世界測地系緯度経度）
　-t_srs EPSG:2448はアウトプットファイルのEPSGコード（大阪なので平面直角6系）
　-tr 10 10 はアウトプットファイルのメッシュサイズ（10m×10m）
　-r bilinear　は投影変換時の内挿方法（バイリニア）
　-dstnodata -9999はnodataな部分（海とか）の値
　union.bil（インプットファイル名）
　union_xy.tif（アウトプットファイル名）

３）GridASCIIにフォーマット変換
そうすると、平面直角6系の投影変換されたGeoTiff形式のDEM（union_xy.tif）ができあがります。残念ながらMapinfo＋vertical mapperではGeiTiff DEMをインポートできないため、Mapinfoでも読めるようにGridASCII（Arc/Info ASCII Grid）に変換してあげます。これもGDALのコマンド（gdal_translate）でできます。

gdal_translate -of AAIGrid union_xy.tif union_xy.asc

という感じ。
本当はgdalwarpの時点で投影変換しながらフォーマットの変換もしちゃったほうがスマートなので、gdalwarpに-of AAIGridオプションをつけてみたのだけれど、「Output driver `AAIGrid' not recognised or does not support direct output file creation.」ということで、gdalwarpで直接GridASCIIを出力するのはできないみたいです。
これで、平面直角6系の10mメッシュがMapinfoにインポートできる形で用意できました。実際にはいちいちGDALのコマンドを手入力するのも大変なのでバッチファイルか何かで実行してあげた方がいいですね。
できあがったunion_xy.ascをMapinfoにインポートすると↓のような感じ。

今回は13個のファイルを処理したので、地図センターに頼むと3.640円！　調べてよかったー。