天気の子、その影響は10メートルの水没⁉︎ その時日本はどうなっていたのか
天気の子の作中で、東京のどの範囲が浸水していたのか。そしてそのときその他の地域、世界はどうなっていたのか。標高データを用いて浸水範囲を考えてみました。
2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html
1. はじめに
2019年7月19日に公開された新海誠監督の最新作映画「天気の子」のBlu-ray&DVDが2020年5月27日に満を持して販売開始！
劇中では、東京の大部分が大雨によって浸水してしまった描写がありました。
もしも、天気の子の舞台であった東京と同程度世界的に海水面が上昇したらどうなったのか……。そんな疑問を持った宙畑編集部は、日本全国、また、世界に広げて衛星データを用いて天気の子の浸水想定範囲を再現してみました。
なお、作中では「降り止まない雨に従来の排水機能が間に合わず、二年以上をかけてゆっくりと海に沈んでいった」との記載があります。本記事で想定するような、「世界的に海水面が上がった状態」に作中ではなっていませんが、もしも東京都と同程度海水面が上昇したら、世界はどうなってしまうのだろうかということを、思考実験的に行った結果になります。
本記事を通して、天気の子は地理的にも的確に描画されていたのだな、と言うことや標高情報面白いから自分でもいろんな地域を細かく見てみよう、などと思って頂けると幸いです。
※本記事は作品のネタばれを一部含みます
※あくまで推測なので、作品の状態を正確に表すものではない、作品の設定に対して意見がある訳ではない、という点を予めご了承ください
2. 本記事で実践する衛星データ活用の概要
天気の子では、物語の終盤に東京の1/3が水没していました。
小説版では、「レインボーブリッジは水に沈み、四本の柱だけが意味ありげな塔のように～」や「新宿駅から山手線に乗った。山手線は今では環状線ではなくなっていて、水没地域を挟んでC字型に途切れている。それぞれの先端にあたる巣鴨駅と五反田駅～」という記述があります。
これらの情報を基に、標高何メートル以下が浸水したのか推測してみます。
最初に、山手線沿線の標高情報を取得し、巣鴨-新宿-五反田沿線は浸水しない閾値を求めます。
標高情報は数値標高モデル、通称DEM(Digital Elevation Model)を利用することで取得できます。衛星や航空機からの観測結果を基に、各地域の標高モデルが生成されています。
詳細に標高情報を取得したいため、まずは国土地理院が公開している航空レーザー測量の結果を用いて標高何メートル以下が浸水したか推定します。
ソースコードは末尾に載せるので、気になる方はご自身の開発環境でも試してみてください。また、国土地理院がブラウザで操作することで標高を取得できるサイトを公開しています。自分の自宅周辺は標高がどれぐらいなのだろうか？と気になった方はぜひご自身でアクセスして確認してみましょう。
23区を中心に、地図と標高情報を重ねた結果を以下に示します。
調整した結果、20m以下の標高の箇所に色付けしています。映画の描写通り、概ね山手線の五反田-巣鴨間の内回りは浸水していることが分かります。
特に、五反田・巣鴨・高田馬場周辺が色づいているため、閾値を決めるべく詳細に確認してみましょう。
まず、五反田周辺の標高情報は以下のようになります。
駅周辺の標高は高いのですが、目黒川沿いは標高が低く、五反田駅自体は10mほどの高さにあるようです。お隣の大崎駅は4mほど。ということで、浸水しないためには、閾値を10m以下に設定する必要があります。
次に、巣鴨駅周辺の標高情報は以下のようになります。
巣鴨駅は20mほどの標高にあることが分かります。お隣の駒込駅は12~15m程度の標高にあることが分かります。
最後に、高田馬場駅周辺の標高情報は以下のようになります。
五反田周辺同様、神田川が流れている影響で標高が低くなっている場所はありますが、山手線沿線は標高が14m以上になっていそうです。
よって、天気の子の描画通りの浸水範囲にするには、五反田駅の条件から10m、巣鴨駅の条件から20m、高田馬場駅周辺の条件から14m以下に閾値を設定すれば良いことが分かります。
試しに10m以下が浸水した、と想定して再度見てみましょう。
標高10m以下が浸水したと仮定すると、五反田と巣鴨は水没せず、映画の描画通りになっていそうです。
ヒロインの陽菜が住んでいる田端駅周辺は、
「右手には緑に茂った桜の木が並び、左手には眼下に広々と眺望が開けている。～そのさらに奥には、雨に濡れた建物がどこまでも続いている。」
という描写があるように、浸水予測範囲との境界に位置しています。
そのため、標高10m以下の地域が浸水したと考えると、概ね作中の設定と同じくなりますが、「レインボーブリッジは水に沈み、四本の柱だけが意味ありげな塔のように～」という記載については、レインボーブリッジ自体の高さが52mほどあるようで、この通りに浸水したとすると山手線沿線すべてが沈んでしまうことになるので、今回はこちらの設定について無視して話を進めます。
※また細かな点では、「東京都の面積の1/3が、今では水の下だった」という描写については、23区の1/3は浸水していそうなことがお分かりいただけるかと思います。
次章では、日本全国の標高10m程度のところまで浸水したとしたらどうなるのか、ということについて調べていきましょう。
3. 日本で標高10m以下の場所をプロットする
標高10mのところまで海面が上昇した場合、日本の各地のどの範囲が浸水することになるのかプロットしていきます。
先ほどは詳細に山手線各駅の標高を知りたかったので、航空レーザーによる測量結果を利用したのですが、今回は広範囲の情報を大まかに知りたいので、衛星による観測結果を用いて描画していきます。今回はTellusに搭載されているASTER GDEM 2というデータを利用します。
早速解析結果ですが、5m以下を青、5mから10mまでの範囲を緑色に示すと以下のようになります。
これだと分かりにくいので、地方ごとに見ていきましょう。
まずは北から。
札幌周辺を中心に、人口が多そうな地域が、標高10m以下となっていそうなことが分かります。実際に、RESASを用いて人口統計情報を取得すると、下記のようになります。上図で色づいている地域と、下記で人口の多い地域に関係がありそうなことが分かりますね。
次に東北地方。
こちらも先ほどと同様に、人が多く住んでいる平野部を中心に、標高が10m以下となっていることが分かります。
人が住んでいる、ということはつまり夜に明るい、ということです。NASAのサイトで確認してみると、なんとなく相関がありそうなことが分かります。
同様に、他の地方も見ていきましょう。海に接している都道府県について、多くの地域の県庁所在地が標高10m以下にありそうだ、ということが分かります。人が住みやすい場所が発展した結果、県庁所在地になっていることを考えると妥当とは思うものの、もう少し内陸部に人口を集中しても良いのかもしれない、などと考えることができますね。
私が気になる地域は詳細にはどうなっているのだろう、と気になった方は、2章を参考にご自身で試してみてください。
4. Tellus以外のツールを使って確認してみる
3章では、Tellusの開発環境を用いて浸水想定範囲を描画してみましたが、Tellus以外のプラットフォームを利用すると、自分でプログラミングせずとも標高、もとい浸水予測範囲を可視化できます。
例えば、Flood Mapというサイトで10m海面が上昇した状態をシミュレーションすると、下記のようになります。
先ほどTellusで得た結果と同じ傾向にあることが分かります。
こちらのサイトでは、日本のみならず、全世界の状態も確認できます。ズームアウトした状態だと変化が分かりにくいですが、ズームインして各国を見ていくと、変化を感じられるかと思います。
余談ですが、国内のサイトだと、産総研が海面上昇シミュレーションシステムというアプリケーションを公開しています。
記事を執筆している5月20日現在では日本域が表示できない状態となっていますが、10m海面が上昇することで世界でどの程度変わるか、お分かり頂けるかと思います。
5.おわりに
天気の子を見て、様々な感想を持たれた方がいるかと思います。
私たち宙畑編集部としては、東京の浸水範囲を見て、「あー、沈んでいるなー」と映画の中のフィクションとして終わらせることなく、異常気象が続いた末に実際にそのような状態に陥ったら私たちの身の回りはどうなってしまうのか、そうならないように身の回りで気を付けられることはなんだろう、などと皆さんが想像する一助となればいいな、と考えて本記事を執筆しました。
近年続く異常気象。温暖化により極域の氷が溶けることで、海水面が上昇している、というニュースをお聞きになった方も多くいらっしゃるかと思います。実際に急に10mも海水面が上昇しなかったとしても、年々海水面は上昇しています。
気象庁のページによると、
世界平均の海面水位は21世紀末（2090～2099年）には、1980～1999年の平均海面水位に対して0.18～0.59m上昇すると予測される。(IPCC第4次評価報告書（2007）)
と記載されています。本記事でご紹介したソースコードやサイトを利用することで、海水面が50cm上昇すると、世界の見た目がどのように変化するのか、試すことができます。
天気の子を見た後に、データとにらめっこしながら、地球の将来について想像し、いま自分にできることは何か、考えてみませんか。
余談
（心の声：宇宙に関する描写もたびたび登場することもあり、新海誠監督の作品が大好きな筆者は、いつか何かしらで衛星データと絡められたらな、と考えていました。劇場で本作を見たときに、これだ！と思い記事化の準備を粛々と進めておりました。）
さて、天気の子のような状態に陥るのは都民である身からするとぜひとも避けたいところですが、作中にも「東京のあの辺はさ、もともとは海だったんだよ。ほんのすこし前ー江戸時代くらいまではね」「江戸そのものがさ、海の入り江だったようだよ。～入り江への戸口が東京だったのさ。～」とあったように、元々はそもそも浸水範囲であった東京の一部地域。
散歩していると水にまつわる史跡も多く見つけられるかと思いますので、ブラタモリじゃありませんが、これを機に、帆高と陽菜の二人が、どこらへんを歩いていたのだろうと想像しながら、昔の地図を片手に歩いてみるのも楽しいかもしれません。
なお、東京港湾事務所のホームページにて、元々はどのようになっていたのか、紹介されているので合わせてぜひご覧ください。
「Tellus」で衛星データを触ってみよう！
日本発のオープン＆フリーなデータプラットフォーム「Tellus」で、まずは衛星データを見て、触ってみませんか？
★Tellusの利用登録はこちらから
ソースコード
2章の都内の標高を詳細に求めるソースコード
こちらのソースコードはTellusの開発環境でなくてもご利用になれます。ご自身の開発環境で実行してみてください。
import requests
import operator
import os
from io import StringIO
from io import BytesIO
from skimage import io
import string
import pandas as pd
import numpy as np
import math
from PIL import Image
import cv2
from colour import Color
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
import matplotlib.pyplot as plt
from matplotlib import cm
今回は国土地理院の標高タイルと標準地図タイルを用います。
地理院タイルについての詳細はこちらのリンクをご覧下さい。
標高タイルにはいくつかの種類がありますが、今回はもっとも精度の高いDEM5Aを使用しました。
地理院タイルについては「Tellusを通してマイ防災マップを作成する」でも触れています。防災マップ作成に興味のある方はそちらの記事も参考にして下さい！
では、まずは標高タイル（DEM5A PNG形式）を取得するための関数を定義します。
def get_hyoko_tile(z, x, y, params={}):
"""
国土地理院標高タイル情報を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
option : dict
APIパスのオプション(entityId)
params : dict
クエリパラメータ
Returns
-------
img: ndarray
タイル画像
http://cyberjapandata.gsi.go.jp/xyz/dem5a/15/29100/12899.txt
"""
url = 'https://cyberjapandata.gsi.go.jp/xyz/dem5a/{}/{}/{}.txt'.format(z, x, y)
r = requests.get(url, params=params)
if r.status_code is not 200:
raise ValueError('status error({}).'.format(r.status_code))
#標高値がない区画は0mに置換する
maptxt = str.replace(r.text, u'e', u'-0.0')
Z = pd.read_csv(StringIO(maptxt), header=None)
Z_array = np.array(Z)
return Z_array
# reference
# http://memomemokun.hateblo.jp/entry/2018/11/01/221545
今回は等高線を用いる代わりに、色分けで標高の変化を捉えます。綺麗なグラデーションを作成するためにはColorを用いるのが楽です。
# 緑とオレンジを呼び出し、それらの色でグラデーションを作成する。
green = Color("green")
orange = Color("orange")
green_orange = list(green.range_to(orange, 14))
green_orange_str = [str(i) for i in green_orange]
作成したリスト、green_orange_strを確認してみます。緑からオレンジに14段階でグラデーションが作成されているのが分かります。
print(green_orange_str)今回は色分けが肝なので、cmapをカスタマイズします。
色分けは18段階（0メートルから18メートル）あり、併せてカラーマップの仕組みに従い、指定した色が0から1の間の数値と対応するようにしています。
cMap = []
for value, colour in zip([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17],\
[green_orange_str[0], green_orange_str[1], green_orange_str[2],\
green_orange_str[3], green_orange_str[4], green_orange_str[5],\
green_orange_str[6], green_orange_str[7], green_orange_str[8],\
green_orange_str[9],green_orange_str[10],green_orange_str[11],\
green_orange_str[12],green_orange_str[13],\
"purple","cyan","red","blue"]):
cMap.append((value/17.0, colour))
customColourMap = LinearSegmentedColormap.from_list("custom", cMap)
作成したもののテストをしてみましょう。ズームレベル、タイルのxy座標を指定し、タイル画像を取得します。
作成した配列情報を描画しますが、この際に先ほど作成したカスタムcmapを用います。
#地形データを読み込む
(z,x,y) = (13,7276,3224) #7279,3229
komagome = get_hyoko_tile(z,x,y)
komagome_ext = np.where((komagome >= 0)&(komagome <= 18),komagome,np.nan)
plt.imshow(komagome_ext,cmap=customColourMap)
plt.colorbar()
plt.show();
画像の確認ができたと思います。指定したタイルによっては
404 Client Error: Not Found for url: http://cyberjapandata.gsi.go.jp/xyz/dem/13/7278/3228.txt
のようなエラーが返ってくることがあります。その部分のタイルは存在しないため、そこを避けるようにタイル画像を取得してください。
今回はベースとなる標準地図タイルと標高タイルを重ね合わせをするため、カラースケール付きのままでは重ね合わせの際に図が大きくずれてしまいます。
そのため、スケールと図を別々に保存する必要があります。以下の手順に従い、カラースケールと地図を分けて保存しましょう。
作成する画像を保存するディレクトリを先に作っておきます。
img_dir_path = "imgs/"
os.makedirs(img_dir_path)
今度は地理院の標準地図タイル画像を取得する関数を定義します。
標準地図の詳細については、こちらをご覧下さい。
def get_base_tile(z, x, y, params={}):
"""
国土地理院標高タイル情報を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
option : dict
APIパスのオプション(entityId)
params : dict
クエリパラメータ
Returns
-------
img: ndarray
タイル画像
https://cyberjapandata.gsi.go.jp/xyz/std/12/3638/1612.txt
"""
url = 'https://cyberjapandata.gsi.go.jp/xyz/std/{}/{}/{}.png'.format(z, x, y)
r = requests.get(url, params=params)
if r.status_code is not 200:
raise ValueError('status error({}).'.format(r.status_code))
return io.imread(BytesIO(r.content))
base = get_base_tile(z,x,y)
plt.figure(figsize=FIGSIZE)
plt.imshow(base)
plt.savefig('imgs/plot_base.png')
# タイル画像の重ね合わせ
background = Image.open("imgs/plot_base.png")
overlay = Image.open("imgs/plot_nobar.png")
background = background.convert("RGBA")
overlay = overlay.convert("RGBA")
new_img = Image.blend(background, overlay, 0.6)
new_img.save("imgs/plot_overlay.png","PNG",dpi=(150,150))
重ね合わせた画像をさらにカラーバーのpng画像と重ね合わせます。二つの画像はサイズがことなるため、そのままでは画像を重ねることができません。
各々の画像の範囲を変えることで、重ね合わせができるようにしましょう。
先ずは画像を読み込み、それぞれの画像サイズを確かめます。
# カラーバーの重ね合わせ
background_new = Image.open("imgs/plot_overlay.png")
overlay_new = Image.open("imgs/plot_onlycbar_tight.png")
background_new = background_new.convert("RGBA")
overlay_new = overlay_new.convert("RGBA")
print(background_new.size)
print(overlay_new.size)
カラーバーをマップの背景画像の右側にくるように配置します。
データのない画像を新しく作成し、そこに背景画像とカラーバー画像を各々配置し、さらにその二つを重ね合わせるという作業をします。
最初に画像の移動を考えて、新しい画像のサイズを決定しましょう。
# 画像の移動
shift_back = (50, 0)
shift_cbar = (432, 82)
# 新しい画像のサイズを計算
nw, nh = map(max, map(operator.add, background_new.size, shift_back), map(operator.add, overlay_new.size, shift_cbar))
print(nw,nh)
準備ができました。上で求めたサイズの画像を新しく用意し、そこに背景画像を貼り付けます。
# new_canvas1にbackground_newを貼り付ける
new_canvas1 = Image.new('RGBA', size=(nw, nh), color=(0, 0, 0, 0))
new_canvas1.paste(background_new,(0,0))
new_canvas1
続けてカラーバーを上で作成した画像に貼り付けてみましょう。
# new_canvas1にoverlay_newを貼り付ける
new_canvas1.paste(overlay_new, (432, 41))
new_canvas1
new_canvas1.save("imgs/plot_overlay_cbar.png","PNG",dpi=(150,150))
これでカラーバーを表示させることができるようになりました。
さらに取得したタイル画像を繋げる関数も定義します。
def get_combined_image(get_image_func, z, topleft_x, topleft_y, size_x=1, size_y=1, params={}):
"""
結合されたタイル画像を取得する
Parameters
----------
get_image_func : function
タイル画像取得メソッド
引数はz, x, y, option, params
z : int
タイルのズーム率
topleft_x : int
左上のタイルのX座標
topleft_y : int
左上のタイルのY座標
size_x : int
タイルを経度方向につなぎわせる枚数
size_y : int
タイルを緯度方向につなぎわせる枚数
option : dict
APIパスのオプション(z,x,y除く)
params : dict
クエリパラメータ
Returns
-------
combined_image: ndarray
結合されたタイル画像
"""
rows = []
blank = np.zeros((256, 256, 3), dtype=np.uint8)
for y in range(size_y):
row = []
for x in range(size_x):
try:
img = get_image_func(z, topleft_x + x, topleft_y + y, params)
except Exception as e:
img = blank
row.append(img)
rows.append(np.hstack(row))
return np.vstack(rows)
ズームレベル、タイル座標を改めて設定し、何枚のタイルを取得するかも指定します。今回は7270から3222を起点にx方向へ8、y方向へ8枚のタイル（合計で2048,2048の配列）を取得します。
# 地理院タイル画像を連続して取得する
(z,x,y) = (13,7270,3222) #7279,3229
size_x = 8
size_y = 8
ベースとなる標準地図を先に描画し、現在の作業ディレクトリにbase.pngとして保存します。
# 描画
plt.figure(figsize=(16,16))
combined_base = get_combined_image(get_base_tile, z, x, y,size_x,size_y)
plt.imshow(combined_base)
plt.savefig('imgs/base.png', dpi=150)
続けて上に重ねる地図を描画し、現在の作業ディレクトリにhyoukou.pngとして保存します。
FIGSIZE = (16,16)
combined_alt = get_combined_image(get_hyoko_tile, z, x, y,size_x,size_y)
combined_alt_ext = np.where((combined_alt >= 0)&(combined_alt <= 10),combined_alt,np.nan)
plt.figure(figsize=FIGSIZE)
mbp = plt.imshow(combined_alt_ext,cmap=customColourMap)
plt.savefig('imgs/hyoukou.png',dpi=150)
# カラーバーのある図を保存（マップの領域は含む）
fig,ax = plt.subplots(figsize=FIGSIZE)
plt.colorbar(mbp,ax=ax)
ax.remove()
# 先ほどの図からバーだけを抜き出した図を保存
plt.savefig('imgs/hyoukou_onlycbar.png',bbox_inches='tight',dpi=150)
二つの画像を重ね、上になるhyoukou.pngは不透明度を下げて表示します。
重ねた画像も同様に保存しておきます。
# タイル画像の重ね合わせ
background = Image.open("imgs/base.png")
overlay = Image.open("imgs/hyoukou.png")
background = background.convert("RGBA")
overlay = overlay.convert("RGBA")
new_img = Image.blend(background, overlay, 0.6)
new_img.save("imgs/overlay_alt.png","PNG",dpi=(150,150))
カラーバーと重ね合わせた画像を結合します。まずは、保存したカラーバーのPNG画像を読み込みましょう。
hyoukou_cbar = Image.open("imgs/hyoukou_onlycbar.png")
hyoukou_cbar = hyoukou_cbar.convert("RGBA")
画像の重ね合わせをするための関数を定義します。
def paste_cbar(background_img,overlay_cbar):
shift_background = (background_img.size[0],0) #背景画像の拡張準備
shift_cbar = (background_img.size[0], background_img.size[1] - background_img.size[1]) #カラーバー画像の拡張準備
# 新しい画像のサイズを計算
img_width, img_height = map(max, map(operator.add, background_img.size, shift_background),\
map(operator.add, overlay_cbar.size, shift_cbar))
# canvas_newにbackground_imgを貼り付ける
canvas_new = Image.new('RGBA', size=(img_width, img_height), color=(0, 0, 0, 0))
canvas_new.paste(background_img,(0,0))
# canvas_newにoverlay_cbarを貼り付ける
canvas_new.paste(overlay_cbar, (background_img.size[0], int((background_img.size[1] - overlay_cbar.size[1])/2)))
return canvas_new #重ね合わせた画像を返す
new_hyoukou_map = paste_cbar(new_img,hyoukou_cbar)
new_hyoukou_map.save("imgs/tokyo_overlay_cbar.png","PNG",dpi=(300,300))
問題なく実行できれば、2章で示したような、標準地図タイルに標高タイル画像が重なった図の隣に、カラーバーが示されているはずです。
次に、各駅を詳細して等高線含めて描画したソースコードを紹介します。
大まかな流れは先ほどまでと同様です。下記を実行することで、巣鴨駅周辺の、等高線が描画された図を取得できます。
# 巣鴨駅〜駒込駅周辺
(z,x,y) = (15,29103,12896)
size_x = 2
size_y = 2
# 描画
plt.figure(figsize=(14,14))
sugamo_base = get_combined_image(get_base_tile, z, x, y,size_x,size_y)
plt.imshow(sugamo_base)
plt.axis('off')
plt.savefig('imgs/sugamo_base.png', dpi=150)
plt.show()
# タイル取得
combined_alt = get_combined_image(get_hyoko_tile, z, x, y,size_x,size_y)
combined_alt_ext = np.where((combined_alt >= 0)&(combined_alt <= 20),combined_alt,np.nan)
#X,Y軸のグリッドを生成
X, Y = np.meshgrid(np.linspace(0,combined_alt_ext.shape[0]-1,combined_alt_ext.shape[0]),\
np.linspace(combined_alt_ext.shape[1]-1,0,combined_alt_ext.shape[1]))
fig = plt.figure(figsize=(14, 14), dpi=150, facecolor='w', edgecolor='k')
#標高に応じて塗りつぶす（0mから20mの範囲で塗りつぶしを行なう）
pcolor = plt.pcolor(X, Y, combined_alt, cmap=customColourMap)
#標高1m間隔で等高線を描く(14〜20m)
elevation = range(14,20,1)
plt.axis('off')
cont = plt.contour(X, Y, combined_alt_ext, levels=elevation, cmap='winter')
cont.clabel(fmt='%1.1f', fontsize=14)
plt.axis('off')
plt.savefig('imgs/sugamo_hyoukou.png', dpi=150)
plt.show()
# タイル画像の重ね合わせ
background = Image.open("imgs/sugamo_base.png")
overlay = Image.open("imgs/sugamo_hyoukou.png")
background = background.convert("RGBA")
overlay = overlay.convert("RGBA")
new_img = Image.blend(background, overlay, 0.6)
new_img.save("imgs/overlay_alt_sugamo.png","PNG",dpi=(150,150))
new_img
3章の全国の標高を大まかに求めるソースコード
TellusではASTER GDEM 2.0の標高データを使うことにより、高さを求めることができます。なお、こちらはTellusの開発環境を用いないと解析できない点を予めご了承ください。
詳しくはASTER GDEM 2.0について解説した下記の記事を参考にしてください。
【ゼロからのTellusの使い方】Jupyter Labで富士山の標高グラフを作成する
まずはライブラリのインポート。そしてトークンとドメインの設定をしましょう。
# ライブラリのインポート
import os, json, requests, math
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
import numpy as np
from progressbar import ProgressBar
import rasterio
from rasterio import features
%matplotlib inline
# トークンとドメインの指定
TOKEN = '自身のトークンを設定してください'
DOMAIN = 'tellusxdp.com'
標高を扱う記事ではおなじみの関数も登場します。この関数により、標高データを求めるためのPNGタイル画像を取得します。
この下から関数の定義が続きますが、実際に引数を操作するものはextract_altのみになります。他の関数は単なるパーツとして考えても問題ありません。
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))
ASTERの標高タイルは、そのままでは単なるRGB値となります。変換式によって標高に戻すための関数を定義しましょう。
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
def rgb_height(image_arr):
"""
一枚の画像に対してcalc_height_chiriin_styleを適用する
（RGB値から標高の計算）
Parameter
---------
image_name: ndarray
"""
pbar = ProgressBar()
heights = []
u = 1
for i in pbar(range(image_arr.shape[0])):
for j in range(image_arr.shape[1]):
R = image_arr[i:i+1,j:j+1,0]
G = image_arr[i:i+1,j:j+1,1]
B = image_arr[i:i+1,j:j+1,2]
heights.append(calc_height_chiriin_style(R, G, B, u))
heights_array = np.array(heights,dtype = "int16")
return heights_array
タイル座標からバウンディングボックス（二つの緯度と二つの経度で作られた領域のこと。通例bboxと略する）を設定する関数を定義します。
def get_bbox(xtile, ytile, zoom):
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)
north_east = num2deg(xtile + 1, ytile, zoom)
south_west = num2deg(xtile, ytile + 1, zoom)
return (south_west[0], south_west[1], north_east[0], north_east[1])
作成したデータをTellus上にpostする関数を定義します。残念なことに一度にpostできるサイズには制限があります。仮に大きなサイズのデータをpostした際にはエラーが生じ、Tellus OS側にはデータはpostされないことになります。
def post_geojson(data):
url = 'https://api.{}/v1/geojson-files'.format(DOMAIN)
headers = {
'Authorization': 'Bearer ' + TOKEN,
'Content-Type': 'application/json'
}
r = requests.post(url, headers=headers, json=data)
while r.status_code is 200:
try:
return r.json()
except ValueError:
print('status error({}).'.format(r.status_code))
次に、使用する色を定義します。色は青から緑へと10段階で変化するようにしています。
上でも書きましたが、postには制限があるため、細かい条件分けを行い一回あたりのpostに含まれるデータ量を減らす必要があります。
postするデータ毎に色を変えれば、条件分けされた高さの判別がしやすくなります。
from colour import Color
blue = Color("blue")
colors = list(blue.range_to(Color("green"),10))
print(colors)
指定した標高を取得する関数を定義します。
zはズームレベル、 x,yはタイル座標、cond_lower、cond_upperは取得する下限と上限になります。
仮に0、10とすれば0メートルより大きく10メートル以下となります。
set_colで色を指定します。デフォルトは青になっています。
def extract_alt(z,x,y,cond_lower,cond_upper,set_col=str(colors[0])):
"""
抽出する高さを指定して、任意の標高タイルから標高を抽出する
その後、抽出した値をポリゴンにし、最終的にgeojsonデータをTellus OS上にpostする
z: zoom level
x: x coordinate
y: y coordinate
cond_lower: the minimum altitudes which you want to retrieve
cond_lower: the maximum altitudes which you want to retrieve
"""
get_alt_tile = get_astergdem2_dsm(z, x, y) #標高タイルを取得
conv_alt_img = rgb_height(get_alt_tile) #標高タイルから高さを求める
conv_alt_img = conv_alt_img.reshape(256,256)
conv_alt_img = np.array(conv_alt_img, dtype = "int16") #ndarrayにし、データタイプは整数値にする
ext_img = (conv_alt_img > cond_lower) & (conv_alt_img <= cond_upper)
shapes = rasterio.features.shapes(conv_alt_img, mask = ext_img) #条件を指定して、画像からシェイプを作成する
# shapesを基にfeaturesを作成
features = list({
'type': 'Feature',
'geometry': s,
'properties': {
'_weight': 4,
'_color': set_col,
'_dashArray': 'null',
'_opacity': 0.5,
'_fillColor': set_col,
'_fillOpacity': 0.2
}
}
for i, (s, v) in enumerate(shapes))
bbox = get_bbox(x, y, z)#ボウンディングボックスを取得
# 作成したfeaturesを地図上に置くために座標値を導く
for feature in features:
coordinates = feature['geometry']['coordinates'][0]
for i in range(len(coordinates)):
coordinates[i] = (coordinates[i][0] * (bbox[2] - bbox[0])/conv_alt_img.shape[0] + bbox[0], \
coordinates[i][1] * (bbox[1] - bbox[3])/conv_alt_img.shape[1] + bbox[3])
# post用のgeojsonを作成
geojson = {
'type': 'FeatureCollection',
'features': features
}
# post用のデータに変換
data = {
'originalName': 'alt_over{}_under_eq{}'.format(cond_lower,cond_upper),
# 'transparency': 0,
'windowNumber': 1,
'data': geojson}
print(post_geojson(data))
では、該当する範囲をポストしていきましょう。
まずは標高が5m以下の範囲をポストします。
(z,x,y) = (6,54,22)
for i in range(4):
for j in range(6):
extract_alt(z,x+i,y+j,0,5,str(colors[0]))
次に、5mから10mまでの範囲をポストします。
(z,x,y) = (6,54,22)
for i in range(4):
for j in range(6):
extract_alt(z,x+i,y+j,5,10,str(colors[5]))
正しく実行できている場合には、3章のようにTellus OSに該当箇所が塗りつぶされて可視化されているかと思います。Tellus OSに描画されていない場合は、画面をリロードしてみてください。
なお、postしたデータをTellus OSから消したい場合は、以下を実行してください。
for i in range(6):
geojsons = get()
print(len(geojsons['items']))
for id in map(lambda item: item['id'], geojsons['items']):
delete(id)
