国土地理院の標高タイルをPythonで取得する

1. この記事について

以前の記事で地理院タイルをJavaScriptで取得する方法を書いた。
国土地理院の標高タイルをJavaScriptで適切に取得する - Qiita
地理院タイルおよび標高タイルの仕様・タイル座標の調べ方などについては上記の記事を参照のこと。

Pythonでも本質的には変わらないが、Pandasを使えば楽なのでここに記す。
また、複数のタイル座標にわたってタイルを取得・結合する方法も紹介する。

もともとは以下の記事で使うために書いた。
PFN製の最適化ツール「Optuna」で富士山を登頂する

2. 標高タイルの取得方法

基本的にはタイル座標を含めたGETリクエストを送るだけ。
Pandasのread_csv関数を使えば、指定したURLからのデータの取得・CSVのパースを一括で行ってくれるので楽になる。
海などで標高にeという値が入っていることがあるが、.replaceメソッドを呼び出すだけで0mに変換できる。

import pandas as pd

def fetch_tile(z, x, y):
    url = "https://cyberjapandata.gsi.go.jp/xyz/dem/{z}/{x}/{y}.txt".format(z=z, x=x, y=y)
    df =  pd.read_csv(url, header=None).replace("e", 0)
    return df.values

返り値の型はnumpy.ndarrayになるようにしている。
これで標高タイルを2次元のNumpy Arrayとして扱うことが可能になる。

3. 使い方

上記のように定義したfetch_tile関数にタイル座標を与えるだけ。
そのままMatplotlibのimshow関数で可視化すると、ちゃんと北が上・東が右になるように表示される(๑˃̵ᴗ˂̵)

nabewari = (13, 7262, 3232) # タイル座標 (z, x, y)
nabewari_tile = fetch_tile(*nabewari)

# 可視化も簡単
import matplotlib.pyplot as plt
plt.imshow(nabewari_tile)

4. 複数タイルの結合

複数のタイル座標にまたがって長方形領域を取得したい場合、以下のような関数を使えば良い。前章のfetch_tileは定義済とする。

def fetch_all_tiles(north_west, south_east):
    """ 北西端・南東端のタイル座標を指定して、長方形領域の標高タイルを取得 """
    assert north_west[0] == south_east[0], "タイル座標のzが一致していません"
    x_range = range(north_west[1], south_east[1]+1)
    y_range = range(north_west[2], south_east[2]+1)
    return  np.concatenate(
        [
            np.concatenate(
                [fetch_tile(north_west[0], x, y) for y in y_range],
                axis=0
            ) for x in x_range
        ],
        axis=1
    )

例えば以下のような領域を取得したい場合

左上と右下のタイル座標を指定すれば良い。
ここではタイル数が2x2だが、長方形領域ならなんでも良い。

tile = fetch_all_tiles((14, 14523, 6464), (14, 14524, 6465))
print(tile.shape) # (512, 512)

以上