はじめに
前回に引き続き、日本発の衛星データプラットフォーム Tellus を Elixir Livebook から使ってみました
前回 までの記事
https://qiita.com/RyoWakabayashi/items/3c78ad939f83a2eb14c4
https://qiita.com/RyoWakabayashi/items/e884ca99f74cd05717ea
今回は降水量データを地図上にオーバーレイしてみます
今回実装したノートブックの全量はこちら
出典
この記事ではTellusの降水観測情報API(試用版)を使用しています
©島津ビジネスシステムズ
実行環境
このリポジトリーの Docker コンテナ上で実行しました
Tellus の使い方
Tellus の使い方は以前の記事を参照してください
降水観測情報の購入（無料版）
Tellus Market から 降水観測情報（無料版） を購入します
トークンの作成
APIで認証するためのトークンを作成します
すでにトークンを作成している場合は同じものを使い回せます
セットアップ
Livebook を起動してノートブックを開きます
以下のコードをノートブック上で実行してください
Mix.install([
{:nx, "~> 0.4"},
{:evision, "~> 0.1"},
{:httpoison, "~> 1.8"},
{:json, "~> 1.4"},
{:kino, "~> 0.8"},
{:kino_maplibre, "~> 0.1.3"}
])
必要なライブラリをインストールします
- nx: 行列演算
- evision: 画像処理
- httpoison: REST API を呼び出す
- json: JSON をエンコード、デコードする
- kino: 実行結果を可視化する
- kino_maplibre: 実行結果を地図上に可視化する
以下のコードを実行するとテキストエリアが表示されるので、 Tellus で作っておいたトークンを入力します
token_input = Kino.Input.text("Token")
TelluSARの商品ID（ Tellus の「購入済み」画面から確認できます）を設定しておきます
# 降水量データの商品ID
rain_product_id = "51f86de1-777f-43e6-821f-a2c1c737cb8b"
降水観測情報API のベースURLを設定しておきます
base_url = "https://sbs.tellus-tools.com"
Tellus の認証
Tellus 用の認証トークン取得処理を定義しておきます
auth_url = "https://www.tellusxdp.com/api/manager/v2/auth/token/"
json_header = {"Content-Type", "application/json"}
get_auth_header = fn product_id ->
request_body = JSON.encode!(%{product_id: product_id})
auth_header = {"Authorization", "Bearer " <> Kino.Input.read(token_input)}
auth_url
|> HTTPoison.post!(request_body, [auth_header, json_header])
|> then(&JSON.decode!(&1.body))
|> then(&{"Authorization", "Bearer " <> &1["token"]})
end
降水量の取得
降水情報 API は緯度経度と日付がパラメーターになっています
get_precipitation = fn lon, lat, ddate ->
auth_header = get_auth_header.(rain_product_id)
"#{base_url}/tellus/wif82af3w39s/rap_af.php?lon=#{lon}&lat=#{lat}&ddate=#{ddate}"
|> HTTPoison.get!([auth_header])
|> then(&JSON.decode!(&1.body))
end
大分市の 2019/9/9 の降水量を取得してみましょう
# 大分市のデータ
precipitation = get_precipitation.(131.62, 33.23, 20190909)
prec の配列が降水量です
配列の長さを確認してみましょう
Enum.count(precipitation["prec"])
配列の長さは 24 、つまり指定した日付の 00:00 〜 23:00 の1時間毎の降水量（mm）を示しています
この API では緯度経度を 0.01 単位で指定できるようになっています
つまり、特定の範囲の降水量を取得する場合、緯度経度を 0.01 ずつずらしながら API を呼び出す必要があります
大分市はあまり雨が降っていないようなので、東京23区のデータを取得してみましょう
東京23区の東端、西端の経度を設定します
# 東京23区の経度範囲
lon_list =
13956..13991
|> Enum.map(&(&1 / 100))
東京23区の北端、南端の経度を設定します
# 東京23区の緯度範囲
lat_list =
3552..3582
|> Enum.map(&(&1 / 100))
何回 API を呼び出すのか見てみましょう
Enum.count(lon_list) * Enum.count(lat_list)
1116 件です
Flow を使って並列処理すると負荷試験のようになってしまうので、 Enum.map で逐次呼び出します
precipitation =
Enum.map(lon_list, fn lon ->
Enum.map(lat_list, fn lat ->
get_precipitation.(lon, lat, 20190909)
|> then(&(&1["prec"]))
end)
end)
それなりに降っていそうです
ただ、このままだとどこにどれだけ降っているのかよく分かりません
可視化してみましょう
降水量の可視化
まず単純に Nx のテンソルに変換してみます
prec_tensor =
precipitation
|> Nx.tensor()
36 * 31 * 24 のテンソルになりました
経度が0.36、緯度が0.31の範囲で24時間分の降水量、ということです
この範囲（東京23区の2019/9/9）での最低降水量は以下のようにして求められます
Nx.reduce_min(prec_tensor)
0 ですね
次に最大降水量を求めてみましょう
テンソルのままだとこの後扱いにくいので、 Float にします
max_prec =
prec_tensor
|> Nx.reduce_max()
|> Nx.to_number()
90.0 mm が最大でした
つまり、この範囲での降水量は 0 〜 90 ということです
地図に重ねて降水量を可視化するとした場合、 0 = 全く降っていないときは透明、 最大の 90 のときは不透明、その中間は半透明にすれば良さそうです
画像は横ピクセル * 縦ピクセル * RGBA（赤緑青透明度） で表すことができます
RGB の部分は雨らしく真っ青として、 A （透明度）を降水量にすると良いですね
ただし、透明度は 0 = 完全な透明から 255 = 完全な不透明の値で表すため 0 - 90 の値をこの幅に合うように調整したほうが見やすくなりそうです
このような操作を正規化と言います
Nx.multiply(tensor, 255 / max_prec) のように、255 をかけて元のテンソルの（最大値 - 最小値）で割れば正規化できます
また、24時間分を一度に表すのは難しいため、まずは 00:00 の降水量だけを見てみます
Nx.slice_along_axis(tensor, 0, 1, axis: 2) のようにすることで、テンソルから特定の範囲（今回であれば特定の時間=axis（軸）のインデックス2について、0からひとつ分）を取り出せます
alpha =
prec_tensor
|> Nx.multiply(255 / max_prec)
|> Nx.slice_along_axis(0, 1, axis: 2)
画像として扱う場合、経度は横ピクセル、緯度は縦ピクセルになります
その幅と高さを取得しておきましょう
{w, h, t} = Nx.shape(prec_tensor)
次に色を真っ青に指定しておきます
今まで RGB と言っておきながらですが、実は Evision (OpenCV) では BGR になっています
なので、真っ青は [255, 0, 0] になります
この色が先ほど取得した幅、高さの分だけ繰り返されるので、 Nx.tile を使って以下のように設定できます
bgr =
[255, 0, 0]
|> Nx.tensor()
|> Nx.tile([w, h, 1])
全面真っ青な画像の各ピクセルについて、透明度を付けます
そのテンソルを Evision のマトリックスにすることで画像として扱えるようになります
更に見やすくするため、縦横10倍にリサイズします
また、大体のところが低い降水量で見にくいため、 Evision.convertScaleAbs でコントラストを強調します
引数の alpha が大きいほどコントラストが強くなります
[bgr, alpha]
|> Nx.concatenate(axis: 2)
|> Evision.Mat.from_nx_2d()
|> Evision.resize({w * 10, h * 10})
|> Evision.convertScaleAbs(alpha: 3)
|> dbg()
降水量は上図のようになりました
青が濃いところはたくさんの雨が降っていて、白いところ（透明なところ）は雨が降っていません
では、これを24時間分表示してみましょう
get_mat = fn hour ->
alpha =
prec_tensor
|> Nx.slice_along_axis(hour, 1, axis: 2)
[bgr, alpha]
|> Nx.concatenate(axis: 2)
|> Evision.Mat.from_nx_2d()
|> Evision.resize({56 * 10, 20 * 10})
|> Evision.convertScaleAbs(alpha: 5)
end
0..23
|> Enum.map(fn t ->
t
|> IO.inspect()
|> get_mat.()
|> Kino.render()
end)
深夜から朝にかけて、西から東に雨雲が通過した様子が見えますね
午後はすっかり晴れていたようです
地図へのオーバーレイ
降水量だけ見ても、どこに降ったのか分かりません
前回同様、地図に重ねてみましょう
地図に画像を表示するためには BASE64 にする必要があります
get_data_url = fn mat ->
Evision.imencode(".png", mat)
|> Base.encode64()
|> then(&"data:image/png;base64,#{&1}")
end
地図表示用の MapLibre にエイリアスを設定します
alias MapLibre, as: Ml
東京23区の範囲を表示する関数を定義します
show_map_overlay = fn img_base64 ->
{left_lon, right_lon, bottom_lat, top_lat} = {139.56, 139.91, 35.52, 35.82}
center = {(left_lon + right_lon) / 2, (bottom_lat + top_lat) / 2}
# タイルの中心を地図の中心にする
Ml.new(
center: center,
zoom: 9,
style: :street
)
# 画像をタイルの座標に配置する
|> Ml.add_source(
"sar-source",
type: :image,
url: img_base64,
coordinates: [
[left_lon, top_lat],
[right_lon, top_lat],
[right_lon, bottom_lat],
[left_lon, bottom_lat],
]
)
# 画像をレイヤーとして地図に重ね、透過する
|> Ml.add_layer(
id: "overlay",
source: "sar-source",
type: :raster,
layout: %{
"visibility" => "visible"
}
)
end
この関数を使って 05:00 の降水量を地図に重ねてみましょう
5
|> get_mat.()
|> get_data_url.()
|> show_map_overlay.()
この時間は川崎市に多くの雨が降ったようです
午後は晴れていたので、午前中だけ各時刻の降水量を地図に重ねます
0..11
|> Enum.map(fn t ->
t
|> IO.inspect()
|> get_mat.()
|> get_data_url.()
|> show_map_overlay.()
|> Kino.render()
end)
しっかり23区のどこにどれだけ降ったのかが分かりますね
まとめ
緯度経度と時刻、測定値を Nx でテンソル化することで、統計や可視化が非常にやりやすくなりました