はじめに
日本発の衛星データプラットフォーム Tellus を使って標高データを取得します
今回は ASTER GDEM 3.0 という全世界の3次元地形データを使用します
ASTER GDEM はバージョン 3.0 から誰でも自由に利用可能になりました
Tellus からも無料でダウンロードできます
いつものように Elixir の Livebook を使っていきます
Livebook は簡単にデータを可視化することができます
実装したノートブックはこちら
出典
Credit: ASTER GDEM is courtesy of METI and NASA
事前準備
Tellus のアカウント作成
Tellus から衛星データを取得するため、まずは Tellus のアカウントを作ります
公式ガイドを参照してください
トークンの作成
Tellus の API を呼び出すためのトークンを作成します
画面右上のログインユーザー名->「開発環境」をクリックします
「トークンの発行」をクリックすると、新しいトークンが作成され、下に表示されます
これで API を呼ぶ準備ができました
実行環境
- Elixir: 1.14.2 OTP 24
- Livebook: 0.8.0
以下のリポジトリーの Docker コンテナ上で起動しています
Docker が使える環境であれば簡単に実行できます
https://docs.docker.com/engine/install/
Docker Desktop を無償利用できない場合は Rancher Desktop を使ってください
私のリポジトリーを使う場合、以下の手順で Livebook を起動できます
git clone https://github.com/RyoWakabayashi/elixir-learning.git
cd elixir-learning
$ docker-compose up
...
Attaching to livebook
livebook | [Livebook] Application running at http://0.0.0.0:<ポート番号>/?token=<認証トークン>
最後に表示される URL をブラウザで開けば Livebook にアクセスできます
右上 New notebook をクリックし、新しいノートブックを開けば準備完了です
セットアップ
必要なモジュールをインストールします
Livebook のセルに以下のコードを入力して実行してください
Mix.install(
[
{:nx, "~> 0.4"},
{:evision, "~> 0.1"},
{:exla, "~> 0.4"},
{:req, "~> 0.3"},
{:kino, "~> 0.8"},
{:kino_maplibre, "~> 0.1"},
{:kino_vega_lite, "~> 0.1"}
],
config: [
nx: [
default_backend: EXLA.Backend
]
]
)
- Nx: 行列演算
- Evision: 画像処理
- EXLA: 行列演算の高速化
- Req: HTTPクライアント
- Kino: Livebook での入出力
- Kino MapLibre: 地図出力
- Kino VegaLite: グラフ出力
また、 Nx による行列演算を EXLA を使って高速化する設定をしています
情報の設定
以下のコードを実行するとテキストエリアが表示されるので、 Tellus で作っておいたトークンを入力します
# Tellus のトークンを入力する
token_input = Kino.Input.password("Token")
Tellus Traveler からデータを探す
データ取得には Tellus Satellite Data Traveler API を使用します
API 仕様は以下を参照
データ取得用モジュール定義
データ取得用のモジュールを定義します
defmodule TellusTraveler do
@base_path "https://www.tellusxdp.com/api/traveler/v1"
@data_path "#{@base_path}/datasets"
defp get_headers(token) do
%{
"Authorization" => "Bearer #{token}",
"Content-Type" => "application/json"
}
end
def get_datasets(token, is_order_required) do
url = "#{@data_path}/?is_order_required=#{is_order_required}"
headers = get_headers(token)
url
|> Req.get!(headers: headers)
|> then(& &1.body["results"])
end
def get_dataset(token, dataset_id) do
url = "#{@data_path}/#{dataset_id}/"
headers = get_headers(token)
url
|> Req.get!(headers: headers)
|> then(& &1.body)
end
def search(token, dataset_id, coordinates) do
url =
if is_list(dataset_id) do
"#{@base_path}/data-search/"
else
"#{@data_path}/#{dataset_id}/data-search/"
end
headers = get_headers(token)
request_body =
%{
intersects: %{
type: "Polygon",
coordinates: coordinates
},
query: %{},
sortby: [
%{
field: "properties.start_datetime",
direction: "asc"
}
]
}
|> Map.merge(
if is_list(dataset_id) do
%{datasets: dataset_id}
else
%{}
end
)
|> Jason.encode!()
url
|> Req.post!(body: request_body, headers: headers)
|> then(& &1.body["features"])
end
def get_data_files(token, dataset_id, data_id) do
url = "#{@data_path}/#{dataset_id}/data/#{data_id}/files/"
headers = get_headers(token)
url
|> Req.get!(headers: headers)
|> then(& &1.body["results"])
end
defp get_data_file_download_url(token, dataset_id, data_id, file_id) do
url = "#{@data_path}/#{dataset_id}/data/#{data_id}/files/#{file_id}/download-url/"
headers = get_headers(token)
url
|> Req.post!(headers: headers)
|> then(& &1.body["download_url"])
end
def download(token, dataset_id, scene_id, dist \\ "/tmp/") do
[dist, scene_id]
|> Path.join()
|> File.mkdir_p()
token
|> get_data_files(dataset_id, scene_id)
|> Enum.map(fn file ->
file_path = Path.join([dist, scene_id, file["name"]])
unless File.exists?(file_path) do
token
|> get_data_file_download_url(dataset_id, scene_id, file["id"])
|> Req.get!(output: file_path)
end
file_path
end)
end
end
データセットの確認
データセットの一覧を取得します
datasets =
token_input
|> Kino.Input.read()
|> TellusTraveler.get_datasets(false)
ダウンロード可能で説明文に「標高」が含まれるデータをテーブルに表示します
elevation_datasets =
datasets
|> Enum.filter(fn dataset ->
dataset["permission"]["allow_network_type"] == "global" &&
String.contains?(dataset["description"], "標高")
end)
Kino.DataTable.new(elevation_datasets, keys: ["name", "id", "description"])
二つのデータセットが取得できました
それぞれ、説明文は以下のようになっています
- ASTER GDEM ver.3
ASTER GDEM 3は、水平解像度30m相当、高さ方向精度約12m(日本域)の標高データです。
経済産業省が開発したASTERセンサのデータをもとに作成されており、値には建物や木の高さを含みます。
一部、データ欠損がある箇所はASTER以外のデータを用いており、建物や木の高さを含まない箇所もあります。
パブリックドメインになりますので著作権表記は不要ですが、クレジットとして「ASTER GDEM is courtesy of METI and NASA」の表記を推奨します。
- AW3D30
AW3D30は、水平解像度30m相当、高さ精度目標5mの標高データで高緯度地域もカバーされております。
こちらのデータはJAXAの「だいち」(陸域観測技術衛星 ALOS)に搭載されたPRISMのデータをもとに作成されています。■ 解像度: 緯度・経度1秒 (30m相当)
■ 高さ精度:目標5m (標準偏差)
■ 関連URL:
- プロダクト概要説明
(https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d_j.htm)- オリジナルデータのフォーマット説明書(https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d30/data/aw3d30v3.2_product_j_e1.2.pdf)※
※TellusOSや公開APIとはフォーマットが異なる場合があります。"
今回は ASTER GDEM の方を使うので、データセット ID に ASTER GDEM の ID を指定します
dataset_id = "3f865d0b-6410-453f-b124-e0bf48544b45"
富士山のシーン選択
富士山のおおよその座標を指定します
mt_fuji_location = {138.73, 35.36}
この座標をぐるっと囲む四角形を用意します
mt_fuji_rectangle = [
[
[elem(mt_fuji_location, 0) - 0.001, elem(mt_fuji_location, 1) - 0.001],
[elem(mt_fuji_location, 0) + 0.001, elem(mt_fuji_location, 1) - 0.001],
[elem(mt_fuji_location, 0) + 0.001, elem(mt_fuji_location, 1) + 0.001],
[elem(mt_fuji_location, 0) - 0.001, elem(mt_fuji_location, 1) + 0.001],
[elem(mt_fuji_location, 0) - 0.001, elem(mt_fuji_location, 1) - 0.001]
]
]
Tellus から、この四角形を含む ASTER GDEM のシーン(ある場所、ある時間のデータ)の一覧を取得します
scenes_list =
token_input
|> Kino.Input.read()
|> TellusTraveler.search(dataset_id, mt_fuji_rectangle)
該当するのは1件だけです
GDEM の仕様を読むと分かりますが、 GDEM は世界中の陸地の座標を緯度経度1度ずつのタイルに分割して保有しています
日時も複数持ってはいないため、狭い範囲を指定すれば基本的には1件だけ該当します
mt_fuji_scene_id =
scenes_list
|> Enum.at(0)
|> then(& &1["id"])
これで mt_fuji_scene_id
は "61f5c40d-b0b6-42ec-8aca-dddeb964b534"
というシーン ID になりました
データのダウンロード
データセット ID とシーン ID を指定してデータをダウンロードします
token_input
|> Kino.Input.read()
|> TellusTraveler.download(dataset_id, mt_fuji_scene_id)
実行結果は以下のように表示されます
["/tmp/61f5c40d-b0b6-42ec-8aca-dddeb964b534/ASTGTMV003_N35E138_dem.tif",
"/tmp/61f5c40d-b0b6-42ec-8aca-dddeb964b534/ASTGTMV003_N35E138_num.tif"]
GDEM の使用によると、 ..._dem.tif
は標高データ、 ..._num.tif
は QA (品質評価)データとなっています
.tif
なので画像です
とりあえず Evision (OpenCV) で開いて表示してみましょう
一点気を付けないといけないのが、通常の画像と違い、色が符号付16bitになっている点です
GDEM は標高データなので、1m単位で表そうとすると、通常の色の範囲 0 〜 255 (符号なし8bit)では足りません
-65535 〜 65535 (符号付16bit) であれば海より低い場所もエベレストのように高い場所もカバーできます
Evision.imread
する際、 flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR()
を指定すれば 16bit 画像も問題なく開くことができます
ただし、表示上 255 より上は全て真っ白になります
"/tmp/#{mt_fuji_scene_id}"
|> File.ls!()
|> Enum.filter(fn filename -> Path.extname(filename) != ".txt" end)
|> Enum.sort()
|> Enum.map(fn filename ->
["/tmp", mt_fuji_scene_id, filename]
|> Path.join()
# 色が 16bit で格納されているため、 IMREAD_ANYDEPTH と IMREAD_ANYCOLOR を指定する
|> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
# 大きすぎるのでリサイズ
|> Evision.resize({640, 640})
end)
|> Kino.Layout.grid(columns: 2)
左下が _dem.tif です
なんとなく右下が低いことは分かりますが、他は真っ白です
このままでは可視化できたとは言えないので分かりやすいように加工しましょう
まず _dem.tif を読み、テンソルにしておきます
mt_fuji_dem =
"/tmp/#{mt_fuji_scene_id}"
|> File.ls!()
|> Enum.find(fn filename -> String.ends_with?(filename, "_dem.tif") end)
|> then(&Path.join(["/tmp", mt_fuji_scene_id, &1]))
|> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
|> Evision.Mat.to_nx(EXLA.Backend)
標高データの最大値と最小値を取り出します
{mt_fuji_min_dig, mt_fuji_max_dig} = {
mt_fuji_dem |> Nx.reduce_min() |> Nx.to_number(),
mt_fuji_dem |> Nx.reduce_max() |> Nx.to_number()
}
結果は以下のようになります
{-5, 3757}
最小値は -5 なので、一部海よりも低い場所があるようです
最大値は 3757 で、これが山頂です
誤差があるため実際の富士山頂とは少し違う値になっています
-5 から 3757 の範囲を 0 から 255 の範囲に変換します
mt_fuji_dem_u8 =
mt_fuji_dem
|> Nx.subtract(mt_fuji_min_dig)
|> Nx.multiply(255 / (mt_fuji_max_dig - mt_fuji_min_dig))
|> Nx.as_type(:u8)
Evision.resize(mt_fuji_dem_u8, {640, 640})
山の稜線が見えました
白い筋が峰、黒い筋が谷です
そして一層白く輝いて見えるのが山頂です
ついでにカラーマップを適用して色を付けてみましょう
COLORMAP_JET
は値の低い順に青 -> 水色 -> 緑 -> 黄 -> 赤と色が変化します
mt_fuji_dem_color = Evision.applyColorMap(mt_fuji_dem_u8, Evision.Constant.cv_COLORMAP_JET())
Evision.resize(mt_fuji_dem_color, {640, 640})
山頂は真っ赤でまるで噴火しているようですね
標高をグラフ化する
続いてグラフにしてみましょう
VegaLite は残念ながら 3D に対応していないため、地球を南北に両断したように、断面で標高を見てみます
断面標高グラフ表示用関数を用意します
display_elevation_at_x = fn dem, x_index, width, max_y ->
plot_data =
dem[[x_index]]
|> Nx.to_flat_list()
|> Enum.with_index()
|> Enum.map(fn {elevation, index} ->
%{
elevation: elevation,
index: index
}
end)
x_scale = %{"domain" => [0, 3601]}
y_scale = %{"domain" => [-500, max_y]}
VegaLite.new(width: width)
|> VegaLite.data_from_values(plot_data)
|> VegaLite.mark(:area)
|> VegaLite.encode_field(:x, "index", type: :quantitative, scale: x_scale)
|> VegaLite.encode_field(:y, "elevation", type: :quantitative, scale: y_scale)
end
画像左端、方位で言えば西端で切ってみましょう
display_elevation_at_x.(mt_fuji_dem, 0, 700, 4000)
こんな感じでグラフは山の形そのものになります
では最も高いところを見てみます
まず最も高い(標高の値が最も大きい)座標を探します
max_index = mt_fuji_dem |> Nx.argmax() |> Nx.to_number()
max_x_index = div(max_index, 3601)
max_y_index = max_index - max_x_index * 3601
{max_x_index, max_y_index}
結果は以下のようになります
{2302, 2618}
この座標の値を確認してみます
mt_fuji_dem[[max_x_index, max_y_index]]
結果は以下のようになります
#Nx.Tensor<
s16
EXLA.Backend<host:0, 0.2220939173.319946764.245111>
3757
>
3757 なので、確かに先ほどの最大値と同じです
ではこの x 座標を指定して富士山をぶった斬りましょう
display_elevation_at_x.(mt_fuji_dem, max_x_index, 700, 4000)
なんと見事に尖っていることでしょう
もちろんグラフ全体の幅によって傾きは変わって見えますが、それにしてもよくこんな山に登る人が大勢いるものだと感心します
グラフをアニメーションにする
一部だけを切り取ってみるのは勿体ないので、 3D ができない代わりにアニメーションで見てみましょう
@GeekMasahiro さんの記事を参考に、グラフを動かします
x_scale = %{"domain" => [0, 3601]}
y_scale = %{"domain" => [-500, 4000]}
widget =
VegaLite.new(width: 700)
|> VegaLite.mark(:area)
|> VegaLite.encode_field(:x, "index", type: :quantitative, scale: x_scale)
|> VegaLite.encode_field(:y, "elevation", type: :quantitative, scale: y_scale)
|> Kino.VegaLite.new()
まずウィジェットを用意すると、空のグラフが表示されます
続いて表示用関数を定義します
animate = fn dem, x_index ->
plot_data =
dem[[x_index]]
|> Nx.to_flat_list()
|> Enum.with_index()
|> Enum.map(fn {elevation, y_index} ->
%{
elevation: elevation,
index: y_index
}
end)
Kino.VegaLite.clear(widget)
Kino.VegaLite.push_many(widget, plot_data)
end
では 0 から 3600 まで値を動かしてみましょう
1 ずつだと凄く時間がかかるので、 20 ずつ動かします
0..3600//20
|> Enum.map(fn x_index ->
animate.(mt_fuji_dem, x_index)
Process.sleep(100)
end)
長いので記事には全編載せませんが(また画像アップロードできなくなるので)、山頂付近では山が隆起したように見えます
富士山と大分市を比較する
では富士山と我が町大分県大分市の標高を比較してみましょう
シーン ID 取得を関数化しておきます
get_scene_id = fn location ->
rectangle = [
[
[elem(location, 0) - 0.001, elem(location, 1) - 0.001],
[elem(location, 0) + 0.001, elem(location, 1) - 0.001],
[elem(location, 0) + 0.001, elem(location, 1) + 0.001],
[elem(location, 0) - 0.001, elem(location, 1) + 0.001],
[elem(location, 0) - 0.001, elem(location, 1) - 0.001]
]
]
scenes_list =
token_input
|> Kino.Input.read()
|> TellusTraveler.search(dataset_id, rectangle)
scenes_list
|> Enum.at(0)
|> then(& &1["id"])
end
大分市のざっくりした座標を指定します
oita_location = {131.64, 33.20}
大分市のシーン ID を取得します
oita_scene_id = get_scene_id.(oita_location)
大分市の GDEM データをダウンロードします
token_input
|> Kino.Input.read()
|> TellusTraveler.download(dataset_id, oita_scene_id)
比較するため、大分市も富士山と同じスケールでカラーマップに変換します
並べて表示してみましょう
oita_dem =
"/tmp/#{oita_scene_id}"
|> File.ls!()
|> Enum.find(fn filename -> String.ends_with?(filename, "_dem.tif") end)
|> then(&Path.join(["/tmp", oita_scene_id, &1]))
|> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
|> Evision.Mat.to_nx(EXLA.Backend)
oita_dem_u8 =
oita_dem
|> Nx.subtract(mt_fuji_min_dig)
|> Nx.multiply(255 / (mt_fuji_max_dig - mt_fuji_min_dig))
|> Nx.as_type(:u8)
oita_dem_color = Evision.applyColorMap(oita_dem_u8, Evision.Constant.cv_COLORMAP_JET())
[mt_fuji_dem_color, oita_dem_color]
|> Kino.Layout.grid(columns: 2)
左下の少し明るいところが九重ですね(これでもかなり高いはずですが)
大分も山だらけではありますが、いかに富士山が突き抜けて高いかが分かります
ではグラフを並べてみましょう
どちらも最も高いところです
oita_max_index = oita_dem |> Nx.argmax() |> Nx.to_number()
oita_max_x_index = div(oita_max_index, 3601)
[
display_elevation_at_x.(mt_fuji_dem, max_x_index, 300, 4000),
display_elevation_at_x.(oita_dem, oita_max_x_index, 300, 4000)
]
|> Kino.Layout.grid(columns: 2)
やはり富士山の尖り具合が強すぎる
こうして見ると地形の違いが直感的に分かりますね
富士山とエベレストを比較する
GDEM には世界中の陸地が入っているので、せっかくなら世界と比べてみましょう
というわけで文字通り世界最高峰のエベレストと比較してみます
mt_everest_location = {86.92, 27.99}
mt_everest_scene_id = get_scene_id.(mt_everest_location)
token_input
|> Kino.Input.read()
|> TellusTraveler.download(dataset_id, mt_everest_scene_id)
富士山よりもエベレストの方がもちろん高いので、今度はエベレストのスケールに他を合わせます
mt_everest_dem =
"/tmp/#{mt_everest_scene_id}"
|> File.ls!()
|> Enum.find(fn filename -> String.ends_with?(filename, "_dem.tif") end)
|> then(&Path.join(["/tmp", mt_everest_scene_id, &1]))
|> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
|> Evision.Mat.to_nx(EXLA.Backend)
{mt_everest_min_dig, mt_everest_max_dig} = {
mt_everest_dem |> Nx.reduce_min() |> Nx.to_number(),
mt_everest_dem |> Nx.reduce_max() |> Nx.to_number()
}
min_dig = Enum.min([mt_fuji_min_dig, mt_everest_min_dig])
max_dig = Enum.max([mt_fuji_max_dig, mt_everest_max_dig])
{min_dig, max_dig}
ヒートマップ化するのは関数にしておきます
get_heatmap = fn dem, min_dig, max_dig ->
dem
|> Nx.subtract(min_dig)
|> Nx.multiply(256 / (max_dig - min_dig))
|> Nx.as_type(:u8)
|> Evision.applyColorMap(Evision.cv_COLORMAP_JET())
end
エベレスト、富士山、大分市を並べます
[
get_heatmap.(mt_everest_dem, min_dig, max_dig),
get_heatmap.(mt_fuji_dem, min_dig, max_dig),
get_heatmap.(oita_dem, min_dig, max_dig)
]
|> Kino.Layout.grid(columns: 3)
エベレストの山頂は右上ですね
山の高さも凄いですが、谷との落差が凄まじいです
エベレスト基準にすると、富士山もかなり大人しくなっています
ましてや大分は全体的に真っ暗、まるで海の中です
次にグラフを並べます
mt_everest_max_x_index =
mt_everest_dem
|> Nx.argmax()
|> Nx.to_number()
|> div(3601)
[
display_elevation_at_x.(mt_everest_dem, mt_everest_max_x_index, 200, 9000),
display_elevation_at_x.(mt_fuji_dem, max_x_index, 200, 9000),
display_elevation_at_x.(oita_dem, oita_max_x_index, 200, 9000)
]
|> Kino.Layout.grid(columns: 3)
エベレストはもう全体的に高すぎます
富士山が丸っと隠れてしまうような圧倒的高さです
しかも富士山以上に先端が尖っています
触れれば確実に刺さります
エベレストと比較してしまうと、大分県内は全域平野のようなものです
海面上昇シミュレーション
標高が分かるということは、例えば海面が現在より10m上がった場合、どこが沈むのかを見ることもできます
@zacky1972 さんが記事にしていますが、 2100 年には 10m 上がっていてもおかしくないのです
大分市で海に沈む場所を可視化してみましょう
標高が 0 (現在の海面) より大きく、 10 以下の場所を抽出します
alpha =
Nx.logical_and(Nx.greater(oita_dem, 0), Nx.less_equal(oita_dem, 10))
|> Nx.select(255, 0)
|> Nx.new_axis(-1)
alpha
|> Nx.as_type(:u8)
|> Evision.Mat.from_nx_2d()
|> Evision.resize({640, 640})
白いところが沈む場所です
海岸線に沿っているので、なんとなく大分県の形が見えますね
上の方に少し見えているのは山口県です
白だと見にくいのと、危機感を出すために赤くして他を透明にします
bgra =
[0, 0, 255]
|> Nx.tensor()
|> Nx.tile([3601, 3601, 1])
|> then(&Nx.concatenate([&1, alpha], axis: 2))
|> Nx.as_type(:u8)
|> Evision.Mat.from_nx_2d()
Evision.resize(bgra, {640, 640})
地図上に表示するため、画像をデータURLに変換する関数を用意します
get_data_url = fn mat ->
Evision.imencode(".png", mat)
|> Base.encode64()
|> then(&"data:image/png;base64,#{&1}")
end
MapLibre で地図上にこの画像を投影します
1 / 3600 上下に足し引きしているのは GDEM の仕様に沿っていますが、緯度方向に少しずれているのは恣意的に調整しています
center = {131.5, 33.5}
coordinates =
[
[131 - 1 / 3600 / 2, 34 - 9 / 3600 / 2],
[132 + 1 / 3600 / 2, 34 - 9 / 3600 / 2],
[132 + 1 / 3600 / 2, 33 - 11 / 3600 / 2],
[131 - 1 / 3600 / 2, 33 - 11 / 3600 / 2],
]
img_base64 = get_data_url.(bgra)
MapLibre.new(center: center, zoom: 8, style: :terrain)
|> MapLibre.add_source("image", type: :image, url: img_base64, coordinates: coordinates)
|> MapLibre.add_layer(id: "overlay", source: "image", type: :raster, paint: %{"raster-opacity" => 0.5})
大分市は海に面しているのもあって、10mも海面が上昇するとかなり沈んでしまいます
そして、残念なことに弊社 oec も海の底です
まとめ
標高データを使うことで、地形を直感的に捉えることができました
また、海面上昇の危機を改めて実感しました
持続可能な社会のために、私たちができることをやっていきましょう
私はとにかく Elixir を推していきます!
@zacky1972 さんの連載記事で、 Elixir で SDGs に貢献できることを示してくれています