1. 이 기사에 대하여

이전 기사에서 지리원 타일을 자바 스크립트로 얻는 방법을 썼다.
자바스크립트에서 국토지리원의 고도 타일을 제대로 얻으세요 - Qiita
지리원 타일 및 고도 타일의 사양 · 타일 좌표의 조사 방법 등에 대해서는 상기의 기사를 참조.

파이썬에서도 본질적으로 변하지 않지만, 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)

이상