能登半島地震の航空写真の一枚画像を作る方法

国土地理院が能登半島地震の航空写真を地理院地図へ公開していた。
地理院地図での可視化例

本記事では、Pythonを使って、地図タイルを一枚のGEOTIFF画像にする方法を記載する。

主に下記記事を参考にした。

必要ライブラリのインポート

import subprocess
import glob
from tqdm import tqdm
from tile_operator.operate import TileOperate, file_to_bounds

範囲の指定方法は、上で引用している記事中にも書いてある通り、https://geojson.io/ を使うと良い(本記事での変数名はfile_path)。

取得する範囲を指定

file_path = "./wajima.geojson"
bbox = file_to_bounds(file_path).bounds()
print(bbox)

to = TileOperate(
    bbox=bbox,
    zoom_level=18,
)

to.set_tile_list()
tile_list = to.get_tile_list()
print(len(tile_list))

地図タイルのURLは下記を参照(本記事での変数名はtile_url)。
https://maps.gsi.go.jp/development/ichiran.html#t20240102noto_wazimahigashi_0102do

地図タイルを取得

tile_name = "20240102noto_wazimanaka_0102do"
tile_url = f"https://cyberjapandata.gsi.go.jp/xyz/{tile_name}" + "/{z}/{x}/{y}.png"

for tile in tqdm(to.tile_list):        
    to.download_tile(tile_url, tile[0], tile[1], output=f"./{tile_name}")
    tile_path = f"./{tile_name}/18/{tile[0]}/{tile[1]}.png"
    to.tile_to_geotiff(tile_path)

取得した地図タイルから一枚の画像を作成

files = []
for y_name in tqdm(glob.glob(f"./{tile_name}/18/*")):
    files = files + glob.glob(f"{y_name}/*.tif")
command = ["gdal_merge.py", "-o", f"./18.tif", "-co", "BIGTIFF=YES"] + files
subprocess.run(command)

作成例

図の上側は輪島市「朝市通り」の周辺で、広範囲に全焼していることが分かる。

本記事が何かのお役に立てば幸いです。