【Python】基盤地図情報のDEMデータをXarrayで読み込んで可視化する

# はじめに
国土地理院が公開している[基盤地図情報](https://www.gsi.go.jp/kiban/)の[数値標高モデル（DEM）](https://fgd.gsi.go.jp/download/ref_dem.html)を[Xarray](https://docs.xarray.dev/en/stable/)で読み込みたいとき、次の2つの手順を踏む必要があります。

1. XMLで提供されるDEMをGeoTIFFに変換する
1. 複数のGeoTIFFファイルを結合し、Xarrayのデータ形式（ `xarray.DataArray` ）として読み込む

手順1は、既に[Pythonを用いた基盤地図情報 数値標高データのgeotiff変換](https://qiita.com/HidKamiya/items/66b1d98503301446460c)という記事などで詳しく取り上げられています。ここでは、主に手順2の方法を紹介します。

# XMLからGeoTIFFへの変換
まず、[基盤地図情報からDEMデータをダウンロード](https://fgd.gsi.go.jp/download/menu.php)し、`PackDLMap.zip` を解凍します。次に、前述の[Pythonを用いた基盤地図情報 数値標高データのgeotiff変換](https://qiita.com/HidKamiya/items/66b1d98503301446460c)で紹介されているPythonスクリプトを用いて、 `PackDLMap.zip` に含まれるXMLデータをGeoTIFFに変換します。ただし、筆者がこの記事のスクリプトをそのまま使用したところ、
```python
AttributeError: module 'numpy' has no attribute 'int'.
```
というエラーが出たので、`np.int` となっている箇所をすべて `int` に書き換えました。

# 複数のGeoTIFFファイルの結合と読み込み
さて、ここからが本記事のメインの内容です。

上のスクリプトを使用すると、`PackDLMap` に含まれるすべての `.zip` を `.tif` に変換することができます。これにより、選択した範囲の2次メッシュごとにGeoTIFFファイルが生成されます。これを結合して `xarray.DataArray` として読み込めるようにします。

## 基盤地図情報のメッシュ
[基盤地図情報のダウンロードページ](https://fgd.gsi.go.jp/download/menu.php)にアクセスすると、次のような画面が表示されます。ここ表示されているメッシュが1次メッシュです。

![FGD_mesh1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2434596/af1ae6ae-eb54-c4f9-a1c7-d4d79a401e34.png)

これを拡大すると下のような画面となり、1次メッシュの中にさらに2次メッシュが表示されます。`PackDLMap.zip` を解凍して得られる XMLデータを含んだ `.zip` ファイルは、この2次メッシュ単位となっています。

![FGD_mesh2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2434596/943140e1-4df6-e95b-d6c8-535b6dcccabf.png)

したがって、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）
- 陸域が含まれないメッシュは欠番になる

すなわち、このメッシュ番号にしたがってファイルをソートし、結合すれば良いことになります。

なお、より正確な命名規則は、[【( ..)φメモメモ】基盤地図情報数値標高モデルのファイル名について](https://qiita.com/jp-96/items/528ff81814b21c6c21e4)を参照してください。基本的にメッシュ番号は緯度・経度に基づいて決まるようなので、上で述べたように南西から北東に向かって増加することになります。

## 結合するためのPythonコード
上の規則にしたがって結合するために、次のような関数を定義します。

```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](https://corteva.github.io/rioxarray/stable/)というライブラリを使用しています。

あとは、`combine_geotiffs_in()` 関数に、GeoTIFFファイルが含まれるディレクトリを指定することで、`xarray.DataArray` が得られます。

## 使用例
具体例として、東京都の10mDEMを描画してみます。

まず、基盤地図情報から、東京都全域が含まれるように2次メッシュを選択します。

![FGD_mesh_tokyo.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2434596/20f4fb97-ca5b-9bf1-9c61-5386beb7da54.png)

この範囲内のDEM10Bのデータをすべてダウンロードし、`PackDLMap.zip` を解凍します。

次に、前述の記事で紹介されているPythonスクリプト（ `zdem2tif.py` ）に修正を施したものを使用して、XMLをGeoTIFFに変換します。

```bash
$ cd PackDLMap
$ python zdem2tif.py *.zip
```

そして、本記事で紹介した関数を使用して、GeoTIFFを結合して読み込みます。

```python
z = combine_geotiffs_in("./PackDLMap")  # パスは適宜変更してください
```

# DEMデータの可視化
では、以上のように `xarray.DataArray` として読み込んだDEMデータを可視化してみます。

## 陸域の地形データ用カラーマップを定義
matplotlibに用意されているカラーマップは優秀ですが、地形の描画に適したものはあまりありません。一応、`terrain` というカラーマップは用意されていますが、平地の細かい地形が同じ色で塗られてしまうという課題があります。そこでまず、陸域の地形データ用のカラーマップを定義します。

```python
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毎に別の色で塗り絵して地図を作ってみる](https://memomemokun.hateblo.jp/entry/2019/02/14/085141)で使用されている色をお借りし、日本域の標高に合わせて陸域用のカラーマップとしました。

## DEMデータの陸域を切り抜く
このままDEMデータを描画すると、標高が0m未満の陸域（海抜ゼロメートル地帯）が海と同じ色で描かれる、あるいは海が標高0mの陸域と同じ色で描かれる、ということになってしまいます。そこで、DEMデータを海岸線が含まれるシェープファイル（ここでは都道府県境界）で切り抜いて陸域のみのデータとします。さらに、陸域の海抜ゼロメートル地帯は一様に標高0mとして扱うことで、もし海に色を塗った場合にも同じ色で塗られないようにします。

```python
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コードの最後の行を次のように変更します。

```python
# 陸域で標高が0m未満の領域は一様に0mとして扱い、海域は元の値を保持する
z_clipped = z_clipped.where(z_clipped >= 0, other=np.where(z_clipped < 0, 0, z))
```

また、都道府県境界のシェープファイルには、[geoBoundaries](https://www.geoboundaries.org/index.html)のデータを使用しました。以下のようにして取得できます。

```bash
$ 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データを描画してみます。投影法には、ランベルト正角円錐図法を使用します。

```python
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()
```

これにより、次のような図が描けます。

![topo_tokyo.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2434596/68d8478b-55f8-0506-56cc-db6e6b323386.png)

# 参考文献
- [Pythonを用いた基盤地図情報 数値標高データのgeotiff変換 - Qiita](https://qiita.com/HidKamiya/items/66b1d98503301446460c)
- [【( ..)φメモメモ】基盤地図情報数値標高モデルのファイル名について - Qiita](https://qiita.com/jp-96/items/528ff81814b21c6c21e4)
- [地球全体を、標高1m毎に別の色で塗り絵して地図を作ってみる - はてなブログ](https://memomemokun.hateblo.jp/entry/2019/02/14/085141)
- [基盤地図情報（数値標高モデル）について - 国土地理院](https://fgd.gsi.go.jp/download/ref_dem.html)