ローカル完結で日本全国の標高を緯度経度指定で取得できる環境を作る
この記事はFOSS4G Advent Calendar 2024 の10日目の記事です
大量の緯度経度に紐づく標高値を取得したい
緯度経度の標高値を取得する場合は外部サービスに頼ると簡単に取得できます。
ですが大量の標高値を取得しようとすると環境構築が大変だったりします。
その環境構築の手順をまとめてみました。
問題
国土地理院さんから緯度経度に紐づく標高値を取得するためのAPIが無料で公開されていますが大量のリクエストに送る事ができません。
公式には以下のように控えるように書かれています。
サーバに過度の負担を与えないでください。
過度の負担を与えると判断したアクセスについて、
国土地理院は予告なく遮断を行う場合があります。
どの程度のリクエストが過度の負荷なのかというと公式の回答では1秒間に10回程度のリクエストが上限とされていて大量のリクエストは控えるように書かれています。
なので、まじめにルールを守って大量の標高値を取得すると時間がかかってしまいます。
最近公式Twitterでこの問題に触れていました。
別の方法でGoogle Elevation APIでも標高値の取得ができますが、1リクエストあたり0.7円かかるので大量のリクエストを送り続けるとお金が水のように溶けてしまいます。
解決方法
ローカルに日本全国の標高モデルを落としこんで緯度経度指定で標高値を取得できる環境を作る。これなら外部サーバーに負荷をかけずにすみます。
標高モデルは国土地理院から無料で公開されているのでそれを使い環境構築を行います。
環境構築
前提条件
OS: Windows 11 Home
WSL: ubuntu 22:04
Anaconda インストール済
PCの空き容量が500GBくらい残っている事
以下のリポジトリをcloneしておいてください。中には環境構築用のヘルパー関数と標高取得用のサンプルコードが含まれています。
# リポジトリのディレクトリ構成
$ tree
├── README.md
├── elevation_env.yml # Anacondaの環境設定ファイル
├── elevation_service.py # 標高取得用のクラス
├── run.py # 実行用のサンプルコード
└── tool
├── cab # 国土地理院からDLしたcabファイルの置き場
│ └── gitignore
├── categorizing.sh # cabファイルを解凍してjgd2000とjgd2011に仕分けるスクリプト ※1
├── xml_jgd2000 # ※1で自動的に仕分けられるxml郡
│ └── gitignore
├── xml_jgd2011 # ※1で自動的に仕分けられるxml郡
│ └── gitignore
├── tif_jgd_to_epsg4326.py # tifのepsgコードを4326に変換するスクリプト
├── merge_tif.py # tifファイルを結合するスクリプト
├── tif_jgd2000 # エコリスのツールで生成したtifの置き場
│ └── gitignore
└── tif_jgd2011 # エコリスのツールで生成したtifの置き場
└── gitignore
Pythonの実行環境を作成
地理空間データの変換にはGdalというライブラリを使用するのですがまじめに手動でインストールするとハマるのでAnacondaの仮想環境を使用します。
# 仮想環境を作成
conda env create -f elevation_env.yml
# gdalがインストールされているかを確認
python3 -c "import osgeo; print('GDAL is installed.')"
標高モデルを取得する
基盤地図情報情報 ダウンロードサービスにアクセスして全国の県単位の標高モデルを落とします。
以下の条件で県単位のダウンロードページまで進めます。
- 検索条件指定
- 10mメッシュ
- 10A（火山基本図の等高線）チェック
- 10B（地形図の等高線）チェック
- 選択方式
- 都道府県または市区町村で選択
- 全国 チェック
- 「選択リストに追加」を押下
- 「ダウンロードファイル確認へ」を押下
| 1 | 2 |
|---|---|
 |
 |
ファイルはcab形式で圧縮されています。
変換前準備
標高モデル内にはjgd2000とjgd2011形式が混ざった状態になっています。jgd2011は東日本大震災の地殻変動を反映した物でjgd2000はそれより前の状態のもので最大6m程度ズレています。
この2つが混ざった状態だとGeoTIFFに変換する際に面倒な事になるので事前に仕分けを行います。
手で仕分けると面倒くさいのでスクリプトで自動的仕分けしてしまいます。
DLしたcabファイル郡をtool/cabディレクトリに配置してスクリプトを実行します。
cd tool
./categorizing.sh
# 仕分けられているか確認。xmlファイルがズラーッと表示される。
ls xml_jgd2000
ls xml_jgd2011
categorizing.sh
# 各種ディレクトリ設定
BASE_DIR=$(pwd) # スクリプト実行ディレクトリ
CAB_DIR="$BASE_DIR/cab" # CABファイルの格納ディレクトリ
DEST_DIR="$BASE_DIR/xml_output" # XMLファイルを移動する先のディレクトリ
TEMP_DIR="$BASE_DIR/temp_extract" # 一時ディレクトリ
# JGDフォルダの作成
JGD2000_DIR="$BASE_DIR/xml_jgd2000"
JGD2011_DIR="$BASE_DIR/xml_jgd2011"
mkdir -p "$CAB_DIR" "$TEMP_DIR" "$DEST_DIR" "$JGD2000_DIR" "$JGD2011_DIR"
# プログレスバーの関数
print_progress_bar() {
local progress=$1
local total=$2
local percent=$((progress * 100 / total))
local bars=$((percent / 2)) # プログレスバーの長さを50にする
local spaces=$((50 - bars)) # 残りのスペース
printf "\r[%-50s] %d%%" "$(printf '#%.0s' $(seq 1 $bars))$(printf ' %.0s' $(seq 1 $spaces))" "$percent"
}
# CABファイルのリスト取得
CAB_FILES=("$CAB_DIR"/*.cab)
TOTAL_CAB=${#CAB_FILES[@]}
CURRENT_CAB=0
# CABファイルの処理
echo "Step 1: Extracting CAB files..."
if [ "$TOTAL_CAB" -eq 0 ]; then
echo "No CAB files found in $CAB_DIR."
else
for CAB_FILE in "${CAB_FILES[@]}"; do
CURRENT_CAB=$((CURRENT_CAB + 1))
print_progress_bar "$CURRENT_CAB" "$TOTAL_CAB"
cabextract -d "$TEMP_DIR" "$CAB_FILE" > /dev/null 2>&1
done
fi
echo ""
# ZIPファイルのリスト取得
ZIP_FILES=($(find "$TEMP_DIR" -name "*.zip"))
TOTAL_ZIP=${#ZIP_FILES[@]}
CURRENT_ZIP=0
# ZIPファイルの処理
echo "Step 2: Extracting ZIP files..."
if [ "$TOTAL_ZIP" -eq 0 ]; then
echo "No ZIP files found in $TEMP_DIR."
else
for ZIP_FILE in "${ZIP_FILES[@]}"; do
CURRENT_ZIP=$((CURRENT_ZIP + 1))
print_progress_bar "$CURRENT_ZIP" "$TOTAL_ZIP"
unzip -d "$TEMP_DIR" "$ZIP_FILE" > /dev/null 2>&1
done
fi
echo ""
# XMLファイルの移動
echo "Step 3: Moving XML files..."
XML_FILES=($(find "$TEMP_DIR" -name "*.xml"))
TOTAL_XML=${#XML_FILES[@]}
CURRENT_XML=0
if [ "$TOTAL_XML" -eq 0 ]; then
echo "No XML files found in $TEMP_DIR."
else
for XML_FILE in "${XML_FILES[@]}"; do
CURRENT_XML=$((CURRENT_XML + 1))
print_progress_bar "$CURRENT_XML" "$TOTAL_XML"
mv "$XML_FILE" "$DEST_DIR"
done
fi
echo ""
# XMLファイルの分別
echo "Step 4: Sorting XML files into JGD2000 and JGD2011..."
cd "$DEST_DIR" || exit
TOTAL_XML_SORT=$(find . -name "*.xml" | wc -l)
CURRENT_XML_SORT=0
if [ "$TOTAL_XML_SORT" -eq 0 ]; then
echo "No XML files to sort."
else
for FILE in *.xml; do
CURRENT_XML_SORT=$((CURRENT_XML_SORT + 1))
print_progress_bar "$CURRENT_XML_SORT" "$TOTAL_XML_SORT"
if grep -q '<gml:Envelope srsName="fguuid:jgd2011.bl">' "$FILE"; then
mv "$FILE" "$JGD2011_DIR"
elif grep -q '<gml:Envelope srsName="fguuid:jgd2000.bl">' "$FILE"; then
mv "$FILE" "$JGD2000_DIR"
fi
done
fi
echo ""
# 一時ディレクトリの削除
echo "Step 5: Cleaning up..."
rm -rf "$TEMP_DIR"
echo "Cleanup complete!"
echo "All steps completed successfully!"
GeoTIFF形式に変換する
pythonから扱いやすくするためにGML形式(xml)からGeoTIFF形式に変換します。
現状ubuntuから変換する方法が存在しないのでwindows用のエコリスさんの標高DEMデータ変換ツールを使い変換します。
このツールはJGD2011とJGD2000が混在していると変換できませんが、事前の仕分け作業を行っているので問題なく変換する事ができます。
windowsとubunutのファイルシステムが違う関係でwindows側にxmlをコピーしてから実行すると高速に変換できるのでtool/xml_jgd2000とtool/xml_jgd2011をwindows側にコピーしてそのディレクトリを指定して変換します。
ツール内の変換結合.vbsを起動して以下のように選択するとtifファイルが作成されます。変換には数時間ほどかかります。
投影法を選択してください。
→ 緯度経度: 0
陰影起状図を作成しますか?
→ いいえ
JPGIS(GML形式)の入っているフォルダを選択してください
→ windows側にコピーしたxml_jgd2000 or xml_jgd2011を選択
海域の標高を選択してください。
→ はい: 0
変換された物の中で結合されたtifとは使用しないのでmerge.tif, mergeLL.vrt, xmlファイル以外のtifファイルをubutnu側のtool/tif_jgd2000と tool/tif_jgd2011ディレクトリにコピーします。
GeoTIFFのEPSGコードを統一する
現状のtifはEPSGコードが4162(JGD2000)のものと6668(JGD2011)のものが混ざっています。
このままだと結合できないので一般的で世界測地系のEPSG:4326に統一します。
cd tool
python3 tif_jgd_to_epsg4326.py
tif_jgd_to_epsg4326.py
import os
from osgeo import gdal
jgd2000_tif_path = './tif_jgd2000'
jgd2011_tif_path = './tif_jgd2011'
def main():
# EPSGコードの定義
jgd2000_epsg = 4612
jgd2011_epsg = 6668
dst_epsg = 4326
# jgd2000 の処理
if os.path.exists(jgd2000_tif_path):
print(f"Processing files in {jgd2000_tif_path}...")
batch_reproject_tiffs(jgd2000_tif_path, jgd2000_epsg, dst_epsg)
# jgd2011 の処理
if os.path.exists(jgd2011_tif_path):
print(f"Processing files in {jgd2011_tif_path}...")
batch_reproject_tiffs(jgd2011_tif_path, jgd2011_epsg, dst_epsg)
def reproject_tiff(input_tiff, output_tiff, src_epsg, dst_epsg):
src_ds = gdal.Open(input_tiff)
if src_ds is None:
print(f"Unable to open {input_tiff}")
return
dst_wkt = f'EPSG:{dst_epsg}'
gdal.Warp(output_tiff, src_ds, dstSRS=dst_wkt, srcSRS=f'EPSG:{src_epsg}')
print(f"Reprojected {input_tiff} to {output_tiff} with EPSG:{dst_epsg}")
def batch_reproject_tiffs(input_dir, src_epsg, dst_epsg):
# 入力ディレクトリ内のすべての TIFF ファイルをループ処理
for filename in os.listdir(input_dir):
if filename.lower().endswith('.tif') or filename.lower().endswith('.tiff'):
input_tiff = os.path.join(input_dir, filename)
output_tiff = os.path.join(input_dir, filename)
reproject_tiff(input_tiff, output_tiff, src_epsg, dst_epsg)
# このファイルが直接実行されたときにのみ main を呼び出す
if __name__ == "__main__":
main()
1つのGeoTIFFに結合する
cd tool
python3 merge_tif.py
merge_tif.py
import os
from osgeo import gdal
def find_tiff_files(directories):
"""複数のディレクトリからTIFFファイルを収集"""
tiff_files = []
for directory in directories:
for root, _, files in os.walk(directory):
for file in files:
if file.endswith('.tif') or file.endswith('.tiff'):
tiff_files.append(os.path.join(root, file))
return tiff_files
def create_vrt(input_files, vrt_file):
"""入力ファイルのリストをVRTに変換"""
gdal.BuildVRT(vrt_file, input_files)
def convert_vrt_to_tiff(vrt_file, output_file):
"""VRTを最終的なTIFFファイルに変換"""
gdal.Translate(output_file, vrt_file)
def main(input_directories, output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
# 複数のディレクトリからTIFFファイルを収集
tiff_files = find_tiff_files(input_directories)
if not tiff_files:
print("No TIFF files found in the specified directories.")
return
vrt_file = os.path.join(output_directory, 'temporary.vrt')
output_file = os.path.join(output_directory, 'merged_all.tif')
# VRT作成とTIFFへの変換
create_vrt(tiff_files, vrt_file)
convert_vrt_to_tiff(vrt_file, output_file)
# 一時的なVRTファイルを削除
os.remove(vrt_file)
print(f"All TIFF files from {len(input_directories)} directories have been merged into {output_file}")
input_directories = [
'./tif_jgd2000',
'./tif_jgd2011'
]
output_directory = './tif_merged'
main(input_directories, output_directory)
結合すると240GBくらいの巨大なtifファイルが作成されます。
結合したtifをQGISで読み込ませて日本の形状が表示されればOKです。
Pythonから緯度経度指定で標高値を取得する
あとは結合したtifをgdalに読み込ませて緯度軽度をするだけです。
標高取得ロジックをラップしたクラスを作成。
from osgeo import gdal, osr
class ElevationService:
def __init__(self, tif_path):
self.dataset = gdal.Open(tif_path, gdal.GA_ReadOnly)
proj = self.dataset.GetProjection()
sr = osr.SpatialReference(wkt=proj)
def get_elevation(self, lat: int, lon: int) -> int | None:
if self.dataset is None:
return None
# 緯度と経度をピクセル座標に変換
gt = self.dataset.GetGeoTransform()
x = int((lon - gt[0]) / gt[1])
y = int((lat - gt[3]) / gt[5])
# ピクセル座標から標高を取得
band = self.dataset.GetRasterBand(1)
elevation = band.ReadAsArray(x, y, 1, 1)[0, 0]
return elevation
def __del__(self):
self.dataset = None
クラスの呼び出し側
import tool.elevation_service as elevation_service
tif_path = "./tool/tif_merged/tif_merged.tif"
elevation_service_ins = elevation_service.ElevationService(tif_path)
lat = 35.6895
lon = 139.6917
elevation = elevation_service_ins.get_elevation(lat, lon)
print(f"lat: {lat}, lon: {lon}, elevation: {elevation}")
緯度経度が取得できればOK.
$ python3 run.py
lat: 35.6895, lon: 139.6917, elevation: 37.900001525878906
おわりに
これで制限を気にする事なく日本中の標高値を取得できる環境ができました。
自分がGIS知識がなく環境構築に時間がかかったので同じような状況にいる方の参考になればうれしいです。
Discussion