はじめに
国土地理院が公開している基盤地図情報の数値標高モデル（DEM）をXarrayで読み込みたいとき、次の2つの手順を踏む必要があります。
- XMLで提供されるDEMをGeoTIFFに変換する
- 複数のGeoTIFFファイルを結合し、Xarrayのデータ形式（
xarray.DataArray）として読み込む
手順1は、既にPythonを用いた基盤地図情報 数値標高データのgeotiff変換という記事などで詳しく取り上げられています。ここでは、主に手順2の方法を紹介します。
XMLからGeoTIFFへの変換
まず、基盤地図情報からDEMデータをダウンロードし、PackDLMap.zip を解凍します。次に、前述のPythonを用いた基盤地図情報 数値標高データのgeotiff変換で紹介されているPythonスクリプトを用いて、 PackDLMap.zip に含まれるXMLデータをGeoTIFFに変換します。ただし、筆者がこの記事のスクリプトをそのまま使用したところ、
AttributeError: module 'numpy' has no attribute 'int'.
というエラーが出たので、np.int となっている箇所をすべて int に書き換えました。
複数のGeoTIFFファイルの結合と読み込み
さて、ここからが本記事のメインの内容です。
上のスクリプトを使用すると、PackDLMap に含まれるすべての .zip を .tif に変換することができます。これにより、選択した範囲の2次メッシュごとにGeoTIFFファイルが生成されます。これを結合して xarray.DataArray として読み込めるようにします。
基盤地図情報のメッシュ
基盤地図情報のダウンロードページにアクセスすると、次のような画面が表示されます。ここ表示されているメッシュが1次メッシュです。
これを拡大すると下のような画面となり、1次メッシュの中にさらに2次メッシュが表示されます。PackDLMap.zip を解凍して得られる XMLデータを含んだ .zip ファイルは、この2次メッシュ単位となっています。
したがって、XMLを変換して得られるGeoTIFFファイルも、この2次メッシュ単位となります。そのファイル名は、例えば FG-GML-5339-45-DEM5A.tif というような形式で、このうち 5339 が1次メッシュの番号、45 が 5339 という1次メッシュ内の2次メッシュの番号を表します。
このメッシュ番号は、上の2枚の図を見ればわかるように、規則的に並んでいます。その規則は、次の通りです。
- 1次メッシュ番号の前2桁（
53）は南から北方向に増加する - 1次メッシュ番号の後2桁（
39）は西から東方向に増加する - 2次メッシュは、1次メッシュを8×8に分割する
- 2次メッシュ番号は、各1次メッシュの南西（左下）の角を基準とする
- 2次メッシュ番号の10の位（
4）は南から北方向に増加する（0〜7） - 2次メッシュ番号の1の位（
5）は西から東方向に増加する（0〜7） - 陸域が含まれないメッシュは欠番になる
すなわち、このメッシュ番号にしたがってファイルをソートし、結合すれば良いことになります。
なお、より正確な命名規則は、【( ..)φメモメモ】基盤地図情報数値標高モデルのファイル名についてを参照してください。基本的にメッシュ番号は緯度・経度に基づいて決まるようなので、上で述べたように南西から北東に向かって増加することになります。
結合するためのPythonコード
上の規則にしたがって結合するために、次のような関数を定義します。
import os
import re
import warnings
from itertools import groupby
import rioxarray
import xarray as xr
FILENAME_PATTERN = re.compile(r"FG-GML-(\d{2})(\d{2})-(\d{2})-DEM\d+[ABC]\.tiff?")
def _parse_filename(filename: str) -> tuple[int, int, int, int]:
match = FILENAME_PATTERN.match(filename)
if not match:
raise ValueError(f"File name is not in the expected format: {filename}")
primarymesh_y, primarymesh_x, secondarymesh = map(int, match.groups())
secondarymesh_y = secondarymesh // 10
secondarymesh_x = secondarymesh % 10
return primarymesh_x, primarymesh_y, secondarymesh_x, secondarymesh_y
def _gather_and_sort_files(directory: str) -> list[tuple[int, int, str]]:
file_info = []
for file in os.listdir(directory):
if FILENAME_PATTERN.match(file):
primarymesh_x, primarymesh_y, secondarymesh_x, secondarymesh_y = (
_parse_filename(file)
)
global_x = primarymesh_x * 8 + secondarymesh_x
global_y = primarymesh_y * 8 + secondarymesh_y
file_info.append((global_x, global_y, file))
return sorted(file_info, key=lambda x: (x[1], x[0]))
def _check_sequentiality(sorted_files: list[tuple[int, int, str]]) -> None:
for i in range(len(sorted_files) - 1):
current = sorted_files[i]
next_file = sorted_files[i + 1]
if not (next_file[0] == current[0] + 1 or next_file[1] == current[1] + 1):
warnings.warn(
f"Specified GeoTIFF files may not be adjacent: {current[2]} and {next_file[2]} "
)
def _combine_files_along_x(
directory: str, file_grouped_by_y: dict[int, list]
) -> dict[int, xr.DataArray]:
combined_along_x = {}
for y, items in file_grouped_by_y.items():
arrays = [
rioxarray.open_rasterio(os.path.join(directory, item[2]), masked=True)
for item in items
]
combined_along_x[y] = xr.concat(arrays, dim="x") # type: ignore
return combined_along_x # type: ignore
def combine_geotiffs_in(directory: str) -> xr.DataArray:
"""
Combine GeoTIFF files converted from XML of Fundamental Geospatial
Data (基盤地図情報) in the specified directory and read them as a single `xarray.DataArray`.
"""
sorted_files = _gather_and_sort_files(directory)
_check_sequentiality(sorted_files)
grouped_by_y = {k: list(v) for k, v in groupby(sorted_files, lambda x: x[1])}
combined_along_x = _combine_files_along_x(directory, grouped_by_y)
combined_along_xy = xr.concat(
[combined_along_x[y] for y in sorted(combined_along_x.keys(), reverse=True)],
dim="y",
)
return combined_along_xy.squeeze()
GeoTIFFを xarray.DataArray として読み込むために、rioxarrayというライブラリを使用しています。
あとは、combine_geotiffs_in() 関数に、GeoTIFFファイルが含まれるディレクトリを指定することで、xarray.DataArray が得られます。
使用例
具体例として、東京都の10mDEMを描画してみます。
まず、基盤地図情報から、東京都全域が含まれるように2次メッシュを選択します。
この範囲内のDEM10Bのデータをすべてダウンロードし、PackDLMap.zip を解凍します。
次に、前述の記事で紹介されているPythonスクリプト（ zdem2tif.py ）に修正を施したものを使用して、XMLをGeoTIFFに変換します。
$ cd PackDLMap
$ python zdem2tif.py *.zip
そして、本記事で紹介した関数を使用して、GeoTIFFを結合して読み込みます。
z = combine_geotiffs_in("./PackDLMap") # パスは適宜変更してください
DEMデータの可視化
では、以上のように xarray.DataArray として読み込んだDEMデータを可視化してみます。
陸域の地形データ用カラーマップを定義
matplotlibに用意されているカラーマップは優秀ですが、地形の描画に適したものはあまりありません。一応、terrain というカラーマップは用意されていますが、平地の細かい地形が同じ色で塗られてしまうという課題があります。そこでまず、陸域の地形データ用のカラーマップを定義します。
import numpy as np
import matplotlib.colors as mcolors
def make_landcmap() -> tuple[mcolors.LinearSegmentedColormap, mcolors.Normalize]:
# 標高とそれに対応する色を定義
elevation_color_pairs = np.array(
[
[0, "#1F4806"],
[60, "#68E36B"],
[95, "#98D685"],
[190, "#F9EFCD"],
[500, "#E0BB7D"],
[600, "#D3A62D"],
[1600, "#997618"],
[1900, "#705B10"],
[2200, "#5F510D"],
[2500, "#A56453"],
[3000, "#5C1D09"],
[3500, "snow"],
]
)
# 標高を0から1の範囲に正規化
elevations = elevation_color_pairs[:, 0].astype(float)
min_elevation = elevations.min()
max_elevation = elevations.max()
normalized_elevations = (elevations - min_elevation) / (
max_elevation - min_elevation
)
# 正規化した標高値と対応する色のペアを作成
normalized_color_pairs = [
[norm, color]
for norm, color in zip(normalized_elevations, elevation_color_pairs[:, 1])
]
# カラーマップを生成
land_cmap = mcolors.LinearSegmentedColormap.from_list(
"land", normalized_color_pairs, N=int(max_elevation - min_elevation + 1)
)
cnorm = mcolors.Normalize(vmin=min_elevation, vmax=max_elevation)
return land_cmap, cnorm
# カラーマップの作成
land_cmap, cnorm = make_landcmap()
なお、このカラーマップの作成にあたっては、地球全体を、標高1m毎に別の色で塗り絵して地図を作ってみるで使用されている色をお借りし、日本域の標高に合わせて陸域用のカラーマップとしました。
DEMデータの陸域を切り抜く
このままDEMデータを描画すると、標高が0m未満の陸域（海抜ゼロメートル地帯）が海と同じ色で描かれる、あるいは海が標高0mの陸域と同じ色で描かれる、ということになってしまいます。そこで、DEMデータを海岸線が含まれるシェープファイル（ここでは都道府県境界）で切り抜いて陸域のみのデータとします。さらに、陸域の海抜ゼロメートル地帯は一様に標高0mとして扱うことで、もし海に色を塗った場合にも同じ色で塗られないようにします。
import geopandas as gpd
# 都道府県境界のシェープファイルを読み込み
pref_shape_path = "path_to_shp" # パスは適宜変更してください
prefs = gpd.read_file(pref_shape_path)
prefs = prefs.to_crs(z.rio.crs)
# DEMデータを海岸線が含まれるシェープファイル（都道府県境界）でクリップ
z_clipped = z.rio.clip(prefs.geometry, prefs.crs)
# 陸域で標高が0m未満の領域は一様に0mとして扱い、海域はNaNとする
z_clipped = z_clipped.where(z_clipped >= 0, other=np.where(z_clipped < 0, 0, np.nan))
今回は、海域をすべて NaN として扱い、描画時に背景色を変更することで海域を表現するようにしています。データに含まれる海域を残したい場合には、上のPythonコードの最後の行を次のように変更します。
# 陸域で標高が0m未満の領域は一様に0mとして扱い、海域は元の値を保持する
z_clipped = z_clipped.where(z_clipped >= 0, other=np.where(z_clipped < 0, 0, z))
また、都道府県境界のシェープファイルには、geoBoundariesのデータを使用しました。以下のようにして取得できます。
$ curl -OL "https://www.github.com/wmgeolab/geoBoundaries/raw/9469f09/releaseData/gbOpen/JPN/ADM1/geoBoundaries-JPN-ADM1-all.zip"
$ unzip -d JPN-prefectures geoBoundaries-JPN-ADM1-all.zip
データの描画
では、DEMデータを描画してみます。投影法には、ランベルト正角円錐図法を使用します。
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
# 地図の中心経度
clon = round(((z_clipped.x.min().values + z_clipped.x.max().values) / 2), 1)
# Figure, Axesオブジェクトの作成
fig, ax = plt.subplots(
subplot_kw={
"projection": ccrs.LambertConformal( # ランベルト正角円錐図法
central_longitude=clon, # 中心経度
standard_parallels=(25, 45), # 標準緯度
)
},
)
ax.set_facecolor('#B7E5FA') # 背景を海域の色にする
# DEMデータの描画
z_clipped.plot(
transform=ccrs.PlateCarree(),
cmap=land_cmap,
norm=cnorm,
cbar_kwargs={
"label": "Elevation (m)",
"orientation": "horizontal",
"pad": 0.18,
"extend": "neither",
},
) # type: ignore
# 都道府県境界を描画（都道府県境界データに海岸線データも含まれている）
ax.add_geometries(
geoms=prefs.geometry,
crs=ccrs.PlateCarree(),
ec="black",
fc="none",
lw=0.5,
)
# タイトルの設定
ax.set_title("Topography of Tokyo (DEM10B)")
# グリッド線の設定
gl = ax.gridlines(
lw=0.5, ls="dashed", color="k", draw_labels=True, x_inline=False, y_inline=False
)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {"rotation": 45}
plt.show()
これにより、次のような図が描けます。
