標高データ（ASTER GDEM 2.0）を取得する

本記事ではJupyterLabで富士山の標高グラフを作成する方法を紹介します。

TellusではASTER GDEM 2.0という標高データを利用することができます。ASTER GDEMとは、NASAの地球観測衛星Terraに搭載されたASTERという光学センサの観測データを元に作成された数値標高モデルです。ASTER GDEM 2.0はそのバージョン2にあたり、緯度経度が1秒角（およそ30メートル）、高さ方向が1メートルの分解能で地球全域の標高を得ることができます。

※標高は建物や木の高さを含んだ数値です。一部、建物や木の高さを含まない箇所もあります。

データはAPIによりPNG標高タイルという形式で取得することができます。PNG標高タイルとは、地図タイルと同一の座標を用いて、PNG画像として標高データを表現するフォーマットです。各ピクセルの緯度経度は左上のピクセル座標を代表値とし、標高はRGB値により求めることが可能です。

APIの使い方は商品情報のAPIリファレンスにまとめられています。

※データを利用する際は利用ポリシーを確認して規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/astergdem2/dsm/{z}/{x}/{y}.png

Name	Input	Description
ｚ	必須	Zoom Level 0-12
ｘ	必須	タイルのX座標
ｙ	必須	タイルのY座標

関数を定義して、APIからASTER GDEM2のPNG標高タイル画像を取得しましょう。

import requests, json
from skimage import io
from io import BytesIO
%matplotlib inline
TOKEN = 'ここには自分のアカウントのトークンを貼り付ける'
def get_astergdem2_dsm(zoom, xtile, ytile):
    """
    ASTER GDEM2 の標高タイルを取得
    Parameters
    ----------
    zoom : str 
        ズーム率 
    xtile : str 
        座標X 
    ytile : number 
        座標Y
    Returns
    -------
    img_array : ndarray
        標高タイル(PNG)
    """
    url = " https://gisapi.tellusxdp.com/astergdem2/dsm/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))
img_12_3626_1617 = get_astergdem2_dsm(12, 3626, 1617)
io.imshow(img_12_3626_1617)

ASTER GDEM提供：METI and NASA

サンプルコードでは富士山の山頂付近の標高タイル（z=12, x=3626, y=161）を取得しています（赤く警告が出る場合がありますが、問題ありません）。

PNG標高タイルデータを利用すれば、三次元の画像として描写することもできます。PNG標高タイル画像は標高をRGB値に変換してありますので、これを標高のデータに戻してみましょう。PNG 標高タイル ―Web 利用に適した標高ファイルフォーマットの考案と実装― PNG Elevation Tile ― を参考しますと、標高は以下の式で求めることが可能です。

x = 2^16*R + 2^8*B + B

(i) x< 2^23の場合

height = x*u

(ii) x = 2^23の場合

height = NA（欠損値）

(iii) x > 2^23の場合

height = (x-2^24)*u

uは標高分解能を表す変数になります。仮に標高分解能が0.01mであれば、-2^23*0.01から2^23*0.01の範囲で標高値を表現できることになります。この式を用いて、RGB値から標高を求める関数を定義しましょう。
＊標高タイルでは海（標高0m）は黒（R=0, G=0, B=0）で表現されます。また、十分な観測ができなかった地点はエラー扱いとなり、赤（R=128, G=0, B=0）で表現されます。

def calc_height_chiriin_style(R, G, B, u=100):
    """
    標高タイルのRGB値から標高を計算する
    """
    hyoko = int(R*256*256 + G * 256 + B)
    if hyoko == 8388608:
        raise ValueError('N/A')
    if hyoko > 8388608:
        hyoko = (hyoko - 16777216)*u
    if hyoko < 8388608:
        hyoko = hyoko*u
    return hyoko
# https://www.jstage.jst.go.jp/article/geoinformatics/26/4/26_155/_pdfを参考にした

関数を利用して標高を計算します。

heights = []
u = 1
for i in range(255):
    for j in range(255):
        R = img_12_3626_1617[i,j,0]
        G = img_12_3626_1617[i,j,1]
        B = img_12_3626_1617[i,j,2]
        try:
            heights.append(calc_height_chiriin_style(R, G, B, u))
        except ValueError as e:
            heights.append(0)
            heights_err.append((i,j))
print(min(heights),max(heights)) #最小標高と最大標高を示す

取得した標高値では最小値が1,329m、最大値が3,753mとなりました。実際の標高(3,776m)よりやや低めの値が求められましたが、近似値としては問題ないかと思われます。

続いて、三次元グラフを描くための準備を行います。

Z_value = np.array(heights)
Z_value2d = Z_value.flatten().reshape(255,255)

それでは実際に描画してみましょう。

import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(0, 255, 1)
y = np.arange(0, 255, 1)
X, Y = np.meshgrid(x, y)
Z = Z_value2d
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x_pos")
ax.set_ylabel("y_pos")
ax.set_zlabel("elevation")
ax.plot_wireframe(X, Y, Z)
plt.show()

ASTER GDEM提供：METI and NASA

上のようなグラフが描画できたかと思います。グラフの見栄えについて様々な設定を行うことも可能です。高さによって色分けをする、グラフの向きを変えるなど、興味のある方は挑戦してください。

富士山の標高を算出する

続いて、富士山の標高をグラフを作ってみましょう。z=12, x=3626, y=1617のタイル、統計138.73度に相当する列に対して、緯度方向のRGB値を読み取り、標高へと変換していきます。

タイル座標からバウンディングボックスを取得します。求める座標の値は左下と左上であり、図にすると以下のようになります。つまり、左下のxy座標はタイル上の最小緯度・経度、右上のxy座標はタイル上の最大緯度・経度となっていることがわかります。

import math
import matplotlib.pyplot as plt
def get_tile_bbox(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """
    タイル座標からバウンディングボックスを取得
    https://tools.ietf.org/html/rfc7946#section-5
    """
    def num2deg(xtile, ytile, zoom):
        # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        return (lon_deg, lat_deg)
    right_top = num2deg(topleft_x + size_x , topleft_y, zoom)
    left_bottom = num2deg(topleft_x, topleft_y + size_y, zoom)
    return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])

関数を実行して、富士山を含んだ地図タイル画像を取得します。

bbox = get_tile_bbox(12, 3626, 1617)
print(bbox)

タイル座標が(138.69140625, 35.31736632923787, 138.779296875, 35.389049966911664)に変換されました。

標高を取得し、視覚化する関数を定義する前に、どのように値が抽出されているのかを確認しておきましょう。

上図のようにdirectionが0であるとき、経度が固定され、緯度方向に標高を計算することになります（上から下へと値を読み取る）。RGB値は各々のピクセル毎に抽出され、標高へと変換されます。directionが1である際にも同様のことを行います。

directionが1では、緯度が固定され、x方向にのみ値が変化します。

標高データは各ピクセルから計算され、それを図示する際にはタイル座標を経度緯度に変換される、というのが以下の関数で行なっていることです。

def plot_height(tile, index, direction, z, x, y, xtile_size=1, ytile_size=1, u=1):
    """
    高度のグラフをプロットする
    Parameters
    ----------
    tile : ndarray
        標高タイル
    bbox : object
        標高タイルのバウンディングボックス
    index : number
        固定するピクセルのインデックス（タイル番号）
    direction : number
        求める標高の方向(0で上から下、1で左から右)
    zoom : number
        ズーム率
    xtile : number
        座標X
    ytile : number
        座標Y
    xtile_size : number
        X方向のタイル数
    ytile_size : number  
        Y方向のタイル数
    u : number
        標高分解能[m]（デフォルト 1.0[m]）
    """
    bbox = get_tile_bbox(z, x, y, xtile_size, ytile_size)
    width = abs(bbox[2] - bbox[0])
    height = abs(bbox[3] - bbox[1])
    if direction == 0:
        line = tile[:,index]
        per_deg = -1 * height/tile.shape[0]
        fixed_deg = width/tile.shape[1] * index + bbox[0]
        deg = bbox[3]
        title = 'height at {} degrees longitude'.format(fixed_deg)
    else:
        line = tile[index,:]
        per_deg = width/tile.shape[1]
        fixed_deg = -1 * height/tile.shape[0] * index + bbox[3]
        deg = bbox[0]
        title = 'height at {} degrees latitude'.format(fixed_deg)
    heights = []
    heights_err = []
    degrees = []
    for k in range(len(line)):
        R = line[k,0]
        G = line[k,1]
        B = line[k,2]
        try:
            heights.append(calc_height_chiriin_style(R, G, B, u))
        except ValueError as e:
            heights.append(0)
            heights_err.append((i,j))
        degrees.append(deg)
        deg += per_deg
    fig, ax= plt.subplots(figsize=(8, 4))
    plt.plot(degrees, heights)
    plt.title(title)
    ax.text(0.01, 0.95, 'Errors: {}'.format(len(heights_err)), transform=ax.transAxes)
    plt.show()

関数を実行し、定めたタイル座標の標高を視覚化します。

plot_height(img_12_3626_1617, 115, 0, 12, 3626, 1617)

上のような図が描けました。山頂の火口付近で200メートルほど標高が下がっている様子がグラフからわかります。

今回のスクリプトです。

import requests, json
from skimage import io
from io import BytesIO
%matplotlib inline
TOKEN = 'ここに自分のトークンを貼り付ける'
def get_astergdem2_dsm(zoom, xtile, ytile):
    """
    ASTER GDEM2 の標高タイルを取得
    Parameters
    ----------
    zoom : str
        ズーム率
    xtile : str
        座標X
    ytile : number
        座標Y
    Returns
    -------
    img_array : ndarray
        標高タイル(PNG)
    """
    url = " https://gisapi.tellusxdp.com/astergdem2/dsm/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))
img_12_3626_1617 = get_astergdem2_dsm(12, 3626, 1617)
io.imshow(img_12_3626_1617)
def calc_height_chiriin_style(R, G, B, u=100):
    """
    標高タイルのRGB値から標高を計算する
    """
    hyoko = int(R*256*256 + G * 256 + B)
    if hyoko == 8388608:
        raise ValueError('N/A')
    if hyoko > 8388608:
        hyoko = (hyoko - 16777216)*u
    if hyoko < 8388608:
        hyoko = hyoko*u
    return hyoko
# https://www.jstage.jst.go.jp/article/geoinformatics/26/4/26_155/_pdfを参考にした
heights = []
u = 1
for i in range(255):
    for j in range(255):
        R = img_12_3626_1617[i,j,0]
        G = img_12_3626_1617[i,j,1]
        B = img_12_3626_1617[i,j,2]
        try:
            heights.append(calc_height_chiriin_style(R, G, B, u))
        except ValueError as e:
            heights.append(0)
            heights_err.append((i,j))
Z_value = np.array(heights)
Z_value2d = Z_value.flatten().reshape(255,255)
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
x = np.arange(0, 255, 1)
y = np.arange(0, 255, 1)
X, Y = np.meshgrid(x, y)
Z = Z_value2d
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x_pos")
ax.set_ylabel("y_pos")
ax.set_zlabel("elevation")
ax.plot_wireframe(X, Y, Z)
plt.show()
def get_tile_bbox(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """
    タイル座標からバウンディングボックスを取得
    https://tools.ietf.org/html/rfc7946#section-5
    """
    def num2deg(xtile, ytile, zoom):
        # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        return (lon_deg, lat_deg)
    right_top = num2deg(topleft_x + size_x , topleft_y, zoom)
    left_bottom = num2deg(topleft_x, topleft_y + size_y, zoom)
    return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])
def plot_height(tile, index, direction, z, x, y, xtile_size=1, ytile_size=1, u=1):
    """
    高度のグラフをプロットする
    Parameters
    ----------
    tile : ndarray
        標高タイル
    bbox : object
        標高タイルのバウンディングボックス
    index : number
        走査するピクセルのインデックス
    direction : number
        走査方向(0で上から下、1で左から右)
    zoom : number
        ズーム率
    xtile : number
        座標X
    ytile : number
        座標Y
    xtile_size : number
        X方向のタイル数
    ytile_size : number  
        Y方向のタイル数
    u : number
        標高分解能[m]（デフォルト 1.0[m]）
    """
    bbox = get_tile_bbox(z, x, y, xtile_size, ytile_size)
    width = abs(bbox[2] - bbox[0])
    height = abs(bbox[3] - bbox[1])
    if direction == 0:
        line = tile[:,index]
        per_deg = -1 * height/tile.shape[0]
        fixed_deg = width/tile.shape[1] * index + bbox[0]
        deg = bbox[3]
        title = 'height at {} degrees longitude'.format(fixed_deg)
    else:
        line = tile[index,:]
        per_deg = width/tile.shape[1]
        fixed_deg = -1 * height/tile.shape[0] * index + bbox[3]
        deg = bbox[0]
        title = 'height at {} degrees latitude'.format(fixed_deg)
    heights = []
    heights_err = []
    degrees = []
    for k in range(len(line)):
        R = line[k,0]
        G = line[k,1]
        B = line[k,2]
        try:
            heights.append(calc_height_chiriin_style(R, G, B, u))
        except ValueError as e:
            heights.append(0)
            heights_err.append((i,j))
        degrees.append(deg)
        deg += per_deg
    fig, ax= plt.subplots(figsize=(8, 4))
    plt.plot(degrees, heights)
    plt.title(title)
    ax.text(0.01, 0.95, 'Errors: {}'.format(len(heights_err)), transform=ax.transAxes)
    plt.show()
plot_height(img_12_3626_1617, 115, 0, 12, 3626, 1617)