標高データってどう使うの？Pythonで地理院の数値標高モデル（DEM）を可視化してみた　～その1：基本編

<gml:tupleList>
地表面,232.30
地表面,233.10
地表面,235.20
...
</gml:tupleList>

<gml:boundedBy>
  <gml:Envelope srsName="fguuid:jgd2011.bl">
    <gml:lowerCorner>34.666666667 136.125</gml:lowerCorner>
    <gml:upperCorner>34.675 136.1375</gml:upperCorner>
  </gml:Envelope>
</gml:boundedBy>

<gml:GridEnvelope>
  <gml:low>0 0</gml:low>
  <gml:high>224 149</gml:high>
</gml:GridEnvelope>

<gml:tupleList>
地表面,232.30
地表面,233.10
地表面,235.20
...
</gml:tupleList>

from lxml import etree
import numpy as np
import rasterio
from rasterio.transform import from_bounds
from pathlib import Path

# --- 入力と出力ファイルを指定 ---
xml_path = Path("FG-GML-XXXXXXXXXXXX.xml")  # GML (XML) ファイルパス（ダウンロードしたファイル名）
output_tif = Path("converted_dem.tif")  # 出力先の GeoTIFF

# --- XMLをパース ---
tree = etree.parse(str(xml_path))
root = tree.getroot()
ns = {"gml": "http://www.opengis.net/gml/3.2"}

# --- グリッドサイズを取得 ---
grid_env = root.find(".//gml:GridEnvelope", ns)
if grid_env is None:
    raise ValueError("GridEnvelope が見つかりません")
low = grid_env.find("gml:low", ns).text.split()
high = grid_env.find("gml:high", ns).text.split()
cols = int(high[0]) - int(low[0]) + 1
rows = int(high[1]) - int(low[1]) + 1

# --- 空間範囲を取得 ---
envelope = root.find(".//gml:Envelope", ns)
lower_corner = list(map(float, envelope.find("gml:lowerCorner", ns).text.split()))
upper_corner = list(map(float, envelope.find("gml:upperCorner", ns).text.split()))

# GMLは lat lon の順なので、ここで逆にして代入
min_lat, min_lon = lower_corner
max_lat, max_lon = upper_corner

# --- 標高データを取得 ---
tuple_list_text = root.find(".//gml:tupleList", ns).text.strip()
elevations_flat = [float(line.split(",")[1]) for line in tuple_list_text.splitlines()]
elevations = np.array(elevations_flat).reshape((rows, cols))

# --- アフィン変換行列を計算 ---
transform = from_bounds(min_lon, min_lat, max_lon, max_lat, cols, rows)

# --- XMLパース処理の直後 ---
print(f"サイズ: {cols} × {rows}")
print(f"範囲: 緯度 {min_lat}〜{max_lat}, 経度 {min_lon}〜{max_lon}")
print(f"標高数: {len(elevations_flat)}")


# --- GeoTIFFとして保存 ---
with rasterio.open(
    output_tif,
    "w",
    driver="GTiff",
    height=rows,
    width=cols,
    count=1,
    dtype=elevations.dtype,
    crs="EPSG:6668",  # JGD2011（日本の測地系）
    transform=transform,
) as dst:
    dst.write(elevations, 1)

print("GeoTIFF変換完了！→", output_tif)


# --- script for test ---
lat = xx.xxxx   # ここに緯度経度（deg）を入れておくと読取りチェックできる
lon = xxx.xxxx

with rasterio.open("converted_dem5a.tif") as src:
    col, row = src.index(lon, lat)  # ← 順番注意
    val = src.read(1)[row, col]
    print(f"標高: {val:.2f} m")

low = grid_env.find("gml:low", ns).text.split()
high = grid_env.find("gml:high", ns).text.split()
cols = int(high[0]) - int(low[0]) + 1
rows = int(high[1]) - int(low[1]) + 1

# --- 標高データを取得 ---
tuple_list_text = root.find(".//gml:tupleList", ns).text.strip()
elevations_flat = [float(line.split(",")[1]) for line in tuple_list_text.splitlines()]
elevations = np.array(elevations_flat).reshape((rows, cols))

transform = from_bounds(min_lon, min_lat, max_lon, max_lat, cols, rows)