はじめに
国土地理院が提供する，「基盤地図情報 数値標高モデル」から指定する地点の標高を読み出す方法についてまとめています．
「その1：基本編」では，ダウンロードした数値標高モデルをGeoTIFFと呼ばれる形式に変換し，指定した緯度・経度の標高を表示するPythonプログラムを紹介しました．
上記に引き続き，「その2：応用編」では，
- 複数のファイル（タイル）を一括で処理してGeoTIFFに変換する
- 使い易いように，緯度・経度を与えて，標高を返す機能を関数化する
- 数値標高モデル（DEM）を，上の画像のように色付き（ヒートマップ）で可視化する
の3つについて紹介したいと思います．プログラムは，ここに貼るには少々長いので，GitHubにここで作ったプログラムをまとめました．ダウンロードして，自由に使ってもらってOKです．
複数のGMLファイルを一括でGeoTIFF形式に変換・インデックス作成
今回扱っている地理院DEMのデータは，GeoTIFFに変換したあともファイル数がかなりの数にのぼります．関西圏だけでも数百枚，それを全国分に拡張しようとすると，数千〜数万ファイルにもなります．そのため，極々狭い範囲を対象にする場合は別ですが，通常はDEMデータを一括で変換するプログラムは必須になると思います．
また，そんな状況で「この座標が含まれるタイルはどれ？」という検索を毎回ファイルを開いて確認していたら，処理がとても重くなってしまいます．そこで，あらかじめ各タイルのカバー範囲（緯度経度の上下左右）を記録したインデックスファイルを作成することにします．
一括で変換，といっても，「その1：基本編」で紹介したプログラムをベースにして，データ格納用ディレクトリ内のデータファイルを順次処理して，変換結果格納用ディレクトリに格納するようにしただけです．
なお，標高データが欠損している場合は，「-9999」をダミーとして挿入する処理を行っています．
- 入力DEMファイル格納ディレクトリ名：input_zips
- 出力TIFFファイル格納ディレクトリ名：output_tifs
output_tifs ディレクトリ内のすべての .tif をスキャンし，緯度経度範囲を辞書化します．
インデックスは，output_tifs/dem_index.jsonにJSON形式で以下のような感じで出力されます．
{
"FG-GML-5236-00-00-DEM5A.tif": {
"min_lon": 136.0,
"min_lat": 34.666666667,
"max_lon": 136.0125,
"max_lat": 34.675
},
"FG-GML-5236-00-01-DEM5A.tif": {
"min_lon": 136.0125,
"min_lat": 34.666666667,
"max_lon": 136.025,
"max_lat": 34.675
},
"FG-GML-5236-00-02-DEM5A.tif": {
"min_lon": 136.025,
"min_lat": 34.666666667,
"max_lon": 136.0375,
"max_lat": 34.675
}
}上の，「comvert_all_dem.py」と，「build_dem_index.py」を連続実行するスクリプトです．以下の手順で，GeoTIFFへの変換とインデックスの作成が一気に終わります．
- 基盤地図情報ダウンロードサービスからダウンロードしたデータファイルを，input_zipsに格納（zipのままでOK）
- pythonで，「dem2tif.py」を実行
地点を指定して標高値を求める，マップで可視化する
緯度，経度を引数として呼び出すと，その地点の標高値を返します．
（スクリプトとして実行すると，プログラム内にべた書きされた座標の標高値を出力します）
対象とする地点の座標が，どのtifファイルに入っているかを，上で作成したインデックスファイルを用いて調べた後，該当するtifファイルから対象とする地点を囲む標高グリッドデータを取得し，双線形補間（Bilinear補間）により標高値を求めます．
リポジトリには「get_elevation.py」というプログラムも入っています．こちらは，バイリニア補間無しで，再近傍のグリッドの標高値をそのまま出力する関数です．
DEMタイルを日本地図上にオーバーレイ表示するプログラムです．
実行すると，「output_tifs」ディレクトリに格納されている.tifファイルのデータをマップ上に可視化します．
実行結果と出力される標高値について
動作テストに使用したデータは，上の，緑で網掛けされた領域のデータで，湖，山地，海のバリエーションがあるものになります．
この領域は，さらに細分化してデータ化されており，
FG-GML-523600-DEM5A-20231113.zip
FG-GML-523601-DEM5A-20231113.zip
FG-GML-523602-DEM5A-20231113.zip
のように
FG-GML-523677-DEM5A-20231113.zip
まで，計74個のzipファイルで提供されています（欠番があったり，DEM5BやCのデータ，AとB双方存在…等々があるのでファイル数とファイル名の番号は一致しません）．
また，例えば1番目の「FG-GML-523600-DEM5A-20231113.zip」には，その領域をさらに細分化して，100個のGML（.xml）ファイルが格納されています．
1つのGMLに対して，1つのTIFファイルが生成されますので，この場合は，計5408個のTIFファイルが生成されていました．
ちなみに，
上記のzipファイルを読み込んで地図上に表示するソフトウェア（基盤地図ビューア）が，国土地理院から配布されています．これを用いると，地図上でマウスクリックすることで，ピンポイントに座標，標高値を確認することができます．
基盤地図ビューアは以下のページから取得できます．
基盤地図情報ダウンロードサービス：https://service.gsi.go.jp/kiban/app
以下では，数点ピックアップして，作成したプログラムから得られる標高と，基盤地図ビューアで表示される標高を比較して，作成したプログラムがちゃんと動いているか検証します．
| 地点 | 座標[deg] | 地理院 基盤地図ビューア | get_elevation （補間無し） | get_elevation_bilinear （補間有り） |
| 1 | Lat:35.19096 Lon:136.11497 | 81.6m | 81.7m | 81.59m |
| 2 | Lat:35.13871 Lon:136.30698 | 730.6m | 729.4m | 730.98m |
| 3 | Lat:35.15118 Lon:136.3915 | 973.1m | 798.5m | 964.8m |
| 4 | Lat:35.17146 Lon:136.44017 | 1142.8m | 1099.8m | 1142.7m |
| 5 | Lat:35.17863 Lon:136.41354 | 1221.2m | 1196.7m | 1222.0m |
標高が低めで平坦な場所から，山地の比較的急峻な場所まで試してみました．平坦な場所や，山地でもそれ程急峻でない（地形図で等高線が詰まっていない）ような場所では，バイリニア補間の有無でそれ程違いはありませんでしたが，急峻な場所ではバイリニア補間有りの方が明らかに地理院の基盤地図ビューアと近い値を求めることができました（赤と青のアンダーライン部分）．
ということで，基盤地図ビューアが必ずしも，実際の地形の標高を正しく表示しているのかどうか，という検証は出来ていませんが，取り敢えず今後しばらくは，バイリニア補間有りのバージョンで色々と処理を実装していきたいと思います．
「plot_dem_overview.py」を用いて表示した可視化プロットは以下のようになりました．
（基盤地図ビューアもそうですが，結構処理が重くて，表示までにかなり時間がかかります）
