ASTERの標高データを取得する
標高データ(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 |
---|---|---|
z | 必須 | Zoom Level 0-12 |
x | 必須 | タイルのX座標 |
y | 必須 | タイルの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)
サンプルコードでは富士山の山頂付近の標高タイル(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()
上のようなグラフが描画できたかと思います。グラフの見栄えについて様々な設定を行うことも可能です。高さによって色分けをする、グラフの向きを変えるなど、興味のある方は挑戦してください。
富士山の標高を算出する
続いて、富士山の標高をグラフを作ってみましょう。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)