前回に引き続き、地図ネタです。
みなさんは「Google Elevation API」をご存知ですか?このAPIを使うと簡単に指定地点の標高が取得できます。さらに2点間を指定した標高データも取得できたりととても便利です。
しかし、Google Maps APIの利用は昨年発表された有料化から少し敬遠している人も多いと思います。(最近制限が緩和されました「Google Geo Developers Blog: Lower pricing and simplified limits with the Google Maps API」)
そこで、Google Maps API以外で標高データを取得する方法はないか調べてみました。
標高データは、国土地理院で国土基盤地図としてデータが公開されています。
※利用法によって申請が必要となる場合がありますので詳細は、国土地理院サイトで確認してください。
こちらのデータを使って何かできないかなと調べていたところ、PostGISラスターがPostGIS2.0から標準で装備されたという記事をみつけましたので早速PostGIS2.0系をダウンロードして使うことにしました。
1.標高データ変換(DEM→GeoTIFF)
国土地理院からダウンロードできる標高データはそのままではPostGISにインポートできないので、一度変換を行う必要がありました。変換ツールを公開しているサイトがありましたので、今回はそちらのツールを使用させていただきました(株式会社エコリスさん、無料公開ありがとうございます。)
上記ツールの使い方はダウンロードして展開後の「変換およびコンパイル方法.txt」を参照して下さい(Windowsのみ対応)。
基本的には、所定フォルダに国土地理院からダウンロードして展開した標高データ(xmlファイル)を格納しておき、変換結合.vbsを実行するだけです。表示されたダイアログで「投影法」と「海域の標高」を聞かれますが、今回はそれぞれ「緯度経度:0」と「はい→0」を選択しました。なお、変換対象の標高データは対象フォルダに入れておけば全て変換対象となりますがその分時間がかかるので注意しましょう。
変換が完了すると「変換終了しました。」というダイアログが表示され、対象フォルダ内にtifファイルが生成されます。
2.GeoTIFF→PostGIS投入SQL生成
次はこの変換したtifファイルからPostGISにデータをインポートするSQLを生成します。これにはPostGISに付属する「raster2pgsql」コマンド使用しました。 コマンドの詳しいオプションは以下に記載されています。
「Using raster2pgsql to load rasters」
で色々と情報をあさって今回はこんなオプションで実行してみました(dem_rasterはテーブル名なので適宜変更してください)。
raster2pgsql -s 4326 -I -M *.tif -F -t 50x50 dem_raster > elev.sql
成功するとカレントディレクトリに「elev.sql」が生成されます。
3.PostGISへデータ投入
PostGISのDBへデータ投入します(ここではDB名をetmplate_postgis_20にしています)。
psql -U postgres -d template_postgis_20 -f elev.sql
データ投入が完了したら、pgAdmin等で該当DBにテーブルが作成されデータが入っていることを確認しましょう。
4. 標高データの取得
さて、いよいよ本題です。指定ポイントの標高データを取得するSQLについて色々調べましたが、明確に記載されているところがなかったので試行錯誤してみたところ、それっぽい値がとれました。
template_postgis_20=# SELECT ST_Value(rast,ST_SETSRID(ST_Point(139.601941,35.885747),4612)) AS Elevation From dem_raster where rast ~ ST_SETSRID(ST_Point(139.60 1941,35.885747),4612); elevation ------------------ 9.80000019073486
5.データ比較
取得できた値が正しいかどうか確認するために、いくつかの地点でGoogle Elevation APIの値と比較してみることにしました。また、最近国土地理院が「標高がわかるWeb地図」を試験公開したのでそちらの値とも比較してみることにしました(今回投入した標高データは10mメッシュのみです)。
場所 | 緯度 | 経度 | 国土地理院 | PostGIS | |
高尾山 | 35.626242 | 139.242056 | 565.0 | 569.0 | 568.0 |
埼玉県庁 | 35.856948 | 139.648847 | 21.8 | 15.8 | 14.9 |
東京都庁 | 35.689521 | 139.691694 | 37.1 | 41.4 | 38.0 |
東京ディズニーランド | 35.635253 | 139.880583 | 8.9 | 5.2 | 7.2 |
東京スカイツリー | 35.709621 | 139.810667 | 6.1 | 1.5 | 0.8 |
各サービスで取得する標高の計算方法の詳細はわかりませんが、取得した値を比較してなんとなく標高取得に成功しているという気がします。あくまで何となくなので実際に何かに利用する場合にはもう少し検討が必要かと思いますが、皆さんも一度試してみてはいかがでしょうか。