Elixir Livebook で衛星データプラットフォーム Tellus から降水量データを取得して地図にオーバーレイする

## はじめに

前回に引き続き、日本発の衛星データプラットフォーム Tellus を Elixir Livebook から使ってみました

https://www.tellusxdp.com/ja/

前回　までの記事

https://qiita.com/RyoWakabayashi/items/3c78ad939f83a2eb14c4
https://qiita.com/RyoWakabayashi/items/e884ca99f74cd05717ea

今回は降水量データを地図上にオーバーレイしてみます

今回実装したノートブックの全量はこちら

https://github.com/RyoWakabayashi/elixir-learning/blob/main/livebooks/tellus/precipitation.livemd

## 出典

この記事ではTellusの降水観測情報API(試用版)を使用しています

©島津ビジネスシステムズ

https://www.shimadzu.co.jp/sbs/

## 実行環境

このリポジトリーの Docker コンテナ上で実行しました

https://github.com/RyoWakabayashi/elixir-learning

## Tellus の使い方

Tellus の使い方は以前の記事を参照してください

https://qiita.com/RyoWakabayashi/items/3c78ad939f83a2eb14c4

## 降水観測情報の購入（無料版）

Tellus Market から 降水観測情報（無料版） を購入します

![スクリーンショット 2022-11-30 23.39.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/5be03847-d6da-adc0-b5a1-2caedc023ceb.png)

## トークンの作成

APIで認証するためのトークンを作成します

すでにトークンを作成している場合は同じものを使い回せます

## セットアップ

Livebook を起動してノートブックを開きます

以下のコードをノートブック上で実行してください

```elixir
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 で作っておいたトークンを入力します

```elixir
token_input = Kino.Input.text("Token")
```

TelluSARの商品ID（ Tellus の「購入済み」画面から確認できます）を設定しておきます

```elixir
# 降水量データの商品ID
rain_product_id = "51f86de1-777f-43e6-821f-a2c1c737cb8b"
```

降水観測情報API のベースURLを設定しておきます

```elixir
base_url = "https://sbs.tellus-tools.com"
```

## Tellus の認証

Tellus 用の認証トークン取得処理を定義しておきます

```elixir
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 は緯度経度と日付がパラメーターになっています

https://www.tellusxdp.com/market/api_reference/51f86de1-777f-43e6-821f-a2c1c737cb8b?lang=ja&for_customer=true

```elixir
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)
```

![スクリーンショット 2022-11-30 23.46.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/4a307d3d-9eb4-d79e-89d9-9ce2a9003583.png)

`prec` の配列が降水量です

配列の長さを確認してみましょう

```elixir
Enum.count(precipitation["prec"])
```

![スクリーンショット 2022-11-30 23.48.22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/5dee7707-90e1-1eea-3d3b-a7e6401855b4.png)

配列の長さは 24 、つまり指定した日付の 00:00 〜 23:00 の1時間毎の降水量（mm）を示しています

この API では緯度経度を 0.01 単位で指定できるようになっています

つまり、特定の範囲の降水量を取得する場合、緯度経度を 0.01 ずつずらしながら API を呼び出す必要があります

大分市はあまり雨が降っていないようなので、東京23区のデータを取得してみましょう

東京23区の東端、西端の経度を設定します

```elixir
# 東京23区の経度範囲
lon_list =
  13956..13991
  |> Enum.map(&(&1 / 100))
```

東京23区の北端、南端の経度を設定します

```elixir
# 東京23区の緯度範囲
lat_list =
  3552..3582
  |> Enum.map(&(&1 / 100))
```

何回 API を呼び出すのか見てみましょう

```elixir
Enum.count(lon_list) * Enum.count(lat_list)
```

![スクリーンショット 2022-11-30 23.56.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/099ed38f-ba9b-71a5-5d64-5f0a827cd969.png)

1116 件です

Flow を使って並列処理すると負荷試験のようになってしまうので、 Enum.map で逐次呼び出します

```elixir
precipitation =
  Enum.map(lon_list, fn lon ->
    Enum.map(lat_list, fn lat ->
      get_precipitation.(lon, lat, 20190909)
      |> then(&(&1["prec"]))
    end)
  end)
```

![スクリーンショット 2022-11-30 23.58.10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/fa99c9e1-337c-1891-0be0-f6d92a4d9378.png)

それなりに降っていそうです

ただ、このままだとどこにどれだけ降っているのかよく分かりません

可視化してみましょう

## 降水量の可視化

まず単純に Nx のテンソルに変換してみます

```elixir
prec_tensor =
  precipitation
  |> Nx.tensor()
```

![スクリーンショット 2022-11-30 23.59.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/823f952b-540c-8a4d-8363-0b4c6f6635aa.png)

36 * 31 * 24 のテンソルになりました

経度が0.36、緯度が0.31の範囲で24時間分の降水量、ということです

この範囲（東京23区の2019/9/9）での最低降水量は以下のようにして求められます

```elixir
Nx.reduce_min(prec_tensor)
```

![スクリーンショット 2022-12-01 0.04.08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b36e4406-417d-adb1-5d1d-159f839a8a73.png)

0 ですね

次に最大降水量を求めてみましょう

テンソルのままだとこの後扱いにくいので、 Float にします

```elixir
max_prec =
  prec_tensor
  |> Nx.reduce_max()
  |> Nx.to_number()
```

![スクリーンショット 2022-12-01 0.05.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c0c3103f-e5d3-3507-09fe-bf513aca5728.png)

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からひとつ分）を取り出せます

```elixir
alpha =
  prec_tensor
  |> Nx.multiply(255 / max_prec)
  |> Nx.slice_along_axis(0, 1, axis: 2)
```

![スクリーンショット 2022-12-01 0.35.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/cdfcf3f5-568f-4625-4117-da72c970c440.png)

画像として扱う場合、経度は横ピクセル、緯度は縦ピクセルになります

その幅と高さを取得しておきましょう

```elixir
{w, h, t} = Nx.shape(prec_tensor)
```

![スクリーンショット 2022-12-01 0.37.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/d380617c-c5e1-697c-3d47-c0c2634b6bf8.png)

次に色を真っ青に指定しておきます

今まで RGB と言っておきながらですが、実は Evision (OpenCV) では BGR になっています

なので、真っ青は `[255, 0, 0]` になります

この色が先ほど取得した幅、高さの分だけ繰り返されるので、 `Nx.tile` を使って以下のように設定できます

```elixir
bgr =
  [255, 0, 0]
  |> Nx.tensor()
  |> Nx.tile([w, h, 1])
```

全面真っ青な画像の各ピクセルについて、透明度を付けます

そのテンソルを Evision のマトリックスにすることで画像として扱えるようになります

更に見やすくするため、縦横10倍にリサイズします

また、大体のところが低い降水量で見にくいため、 `Evision.convertScaleAbs` でコントラストを強調します

引数の `alpha` が大きいほどコントラストが強くなります

```elixir
[bgr, alpha]
|> Nx.concatenate(axis: 2)
|> Evision.Mat.from_nx_2d()
|> Evision.resize({w * 10, h * 10})
|> Evision.convertScaleAbs(alpha: 3)
|> dbg()
```

![スクリーンショット 2022-12-01 0.44.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/df70baf6-ffe4-c181-333a-9b0aad88dd6b.png)

降水量は上図のようになりました

青が濃いところはたくさんの雨が降っていて、白いところ（透明なところ）は雨が降っていません

では、これを24時間分表示してみましょう

```elixir
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
```

```elixir
0..23
|> Enum.map(fn t ->
  t
  |> IO.inspect()
  |> get_mat.()
  |> Kino.render()
end)
```

![スクリーンショット 2022-12-01 0.48.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/aa27f4c9-e4f9-bf04-b3ca-84aa6fb03d27.png)
![スクリーンショット 2022-12-01 0.48.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/2be5ffee-1dfb-4b68-27a5-8b0735456708.png)
![スクリーンショット 2022-12-01 0.49.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c3b7b52e-74ac-0df7-18f5-aff4863c55f2.png)
![スクリーンショット 2022-12-01 0.49.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8936167b-4d9b-6384-b775-c431c3f6dabe.png)
![スクリーンショット 2022-12-01 0.49.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/5e3ae898-c991-87a2-ac44-b5cfc156741f.png)

深夜から朝にかけて、西から東に雨雲が通過した様子が見えますね

午後はすっかり晴れていたようです

## 地図へのオーバーレイ

降水量だけ見ても、どこに降ったのか分かりません

前回同様、地図に重ねてみましょう

地図に画像を表示するためには BASE64 にする必要があります

```elixir
get_data_url = fn mat ->
  Evision.imencode(".png", mat)
  |> Base.encode64()
  |> then(&"data:image/png;base64,#{&1}")
end
```

地図表示用の MapLibre にエイリアスを設定します

```elixir
alias MapLibre, as: Ml
```

東京23区の範囲を表示する関数を定義します

```elixir
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 の降水量を地図に重ねてみましょう

```elixir
5
|> get_mat.()
|> get_data_url.()
|> show_map_overlay.()
```

![スクリーンショット 2022-12-01 0.55.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/89e56699-6055-4c32-6e09-17d178a89d00.png)

この時間は川崎市に多くの雨が降ったようです

午後は晴れていたので、午前中だけ各時刻の降水量を地図に重ねます

```elixir
0..11
|> Enum.map(fn t ->
  t
  |> IO.inspect()
  |> get_mat.()
  |> get_data_url.()
  |> show_map_overlay.()
  |> Kino.render()
end)
```

![スクリーンショット 2022-12-01 0.57.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b8e7ccfd-2830-586e-27dc-7e040db1437b.png)
![スクリーンショット 2022-12-01 0.58.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/3223cd64-3037-165d-07fd-366bf0217f94.png)
![スクリーンショット 2022-12-01 0.58.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/cb0ef7ef-9770-8b93-813d-0a9b40a1902b.png)
![スクリーンショット 2022-12-01 0.59.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/28b66ccf-4bfa-aae7-2dd6-5d54f539132f.png)

しっかり23区のどこにどれだけ降ったのかが分かりますね

## まとめ

緯度経度と時刻、測定値を Nx でテンソル化することで、統計や可視化が非常にやりやすくなりました