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

## はじめに

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

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

前回の記事

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

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

SAR とは何かについては以下を参照してください

https://sorabatake.jp/3364/

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

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

この記事ではASNARO-2のデータを使用しています

Original data provided by NEC

## 参考記事

https://sorabatake.jp/7325/

## 実行環境

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

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

## Tellus の使い方

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

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

## TelluSAR　の購入（無料）

Tellus Market から TelluSAR(API) を購入します

![スクリーンショット 2022-11-15 16.16.37.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/56e432b1-ae1a-cfac-915a-05b3f8e12860.png)

## トークンの作成

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

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

## セットアップ

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

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

```elixir
Mix.install([
  {:evision, "~> 0.1"},
  {:httpoison, "~> 1.8"},
  {:json, "~> 1.4"},
  {:kino, "~> 0.7"},
  {:kino_maplibre, "~> 0.1.3"}
])
```

必要なライブラリをインストールします

- evision: 画像処理
- httpoison: REST API を呼び出す
- json: JSON をエンコード、デコードする
- kino: 実行結果を可視化する
- kino_maplibre: 実行結果を地図上に可視化する

以下のコードを実行するとテキストエリアが表示されるので、 Tellus で作っておいたトークンを入力します

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

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

```elixir
sar_product_id = "6a2ae4fc-2a62-4483-861a-906c716add07"
```

TelluSAR のベースURLを設定しておきます

```elixir
base_url = "https://tellusar.tellusxdp.com"
```

SARデータ取得元の衛星IDを設定します

今回は NEC 所有の ASNARO-2 （あすなろつー）のデータを使用します

```elixir
satellite_id = "asnaro2"
```

https://jpn.nec.com/ad/cosmos/asnaro2/index.html

## Tellus の認証

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

認証トークンの有効期限が短いため、以降 API を呼び出す度に認証トークンを取得し直します

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

## シーン一覧の取得

以下のコードを実行すると、取得できるシーン（日時、座標の組み合わせ）の一覧が取得できます

```elixir
auth_header = get_auth_header.(sar_product_id)

scenes_url = "#{base_url}/api/v2/#{satellite_id}/search"

scenes =
  scenes_url
  |> HTTPoison.get!([auth_header, json_header])
  |> then(&JSON.decode!(&1.body))
```

![スクリーンショット 2022-11-15 16.25.26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/1b2a57ce-53de-559a-1711-c200ef0a1980.png)

取得したシーンの中から東京近郊のものを取り出してみます

```elixir
# 東京近郊のデータ
scenes["data"]["scenes"]
|> Enum.filter(&(&1["left_bottom_lat"] >= 35 && &1["left_bottom_lat"] <= 36))
```

![スクリーンショット 2022-11-15 16.27.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/ebae978c-d613-892c-c2bb-daf6777dcd9b.png)

## シーン情報の取得

取得できたシーン一覧の中から、使いたいシーンを指定します

```elixir
scene_id = "AS200534028728-190103"
```

指定したシーンの情報を取得します

```elixir
auth_header = get_auth_header.(sar_product_id)

scenes_info_url = "#{base_url}/api/v2/#{satellite_id}/search/#{scene_id}"

scenes_info =
  scenes_info_url
  |> HTTPoison.get!([auth_header, json_header])
  |> then(&JSON.decode!(&1.body))
```

![スクリーンショット 2022-11-15 17.09.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/435516c0-e314-45e1-b4b0-546ff283e969.png)

　差分干渉のペアとして使えるシーン（座標がほぼ同じシーン）を取得します

```elixir
auth_header = get_auth_header.(sar_product_id)

scenes_after_url = "#{base_url}/api/v2/#{satellite_id}/search/#{scene_id}/afters"

scenes_after =
  scenes_after_url
  |> HTTPoison.get!([auth_header, json_header])
  |> then(&JSON.decode!(&1.body))
```

![スクリーンショット 2022-11-15 17.10.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/d4c5eadd-13eb-2810-3644-da1d7019d645.png)

## 差分干渉データ作成依頼

差分干渉用のシーンIDを指定します

```elixir
after_scene_id = "AS200555328728-190117"
```

Tellus 上で差分干渉データを作成するように要求します

```elixir
auth_header = get_auth_header.(sar_product_id)

work_url = "#{base_url}/api/v2/works"

work_body = %{
  "satellite" => satellite_id,
  "before_scene_id" => scene_id,
  "after_scene_id" => after_scene_id,
  "polarisation" => "HH",
  "nlook_rg" => 5,
  "nlook_az" => 7,
  "filter" => 0,
  "beam_number" => 1
}

work =
  work_url
  |> HTTPoison.post!(JSON.encode!(work_body), [auth_header, json_header])
  |> then(&JSON.decode!(&1.body))
```

![スクリーンショット 2022-11-15 16.33.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7d7803f5-20c9-b970-6824-a2c3823ae4d4.png)

まだ差分干渉データを作成したことのないペアであれば、 `exist_flag` が false になり、 Tellus 上で処理が開始されます

後で使用するため、レスポンスに含まれる `work_id` を取得しておきます

```elixir
work_id = work["data"]["work_id"]
```

処理状況を確認します

```elixir
auth_header = get_auth_header.(sar_product_id)

work_status_url = "#{base_url}/api/v2/works/#{work_id}"

work_status =
  work_status_url
  |> HTTPoison.get!([auth_header, json_header])
  |> then(&JSON.decode!(&1.body))
```

まだ処理中の場合、 `complete_date` が `nil` になっています

![スクリーンショット 2022-11-15 16.35.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/37bf0a25-851f-7d9d-cf6c-39c35c15d67e.png)

数分後にまた処理状況を確認し、 `complete_date` に値が入っていれば差分干渉データがダウンロードできます

## PNGの取得

SAR データを PNG 形式の画像として取得します

画像の拡大率を指定します

```elixir
z = 11
```

この拡大率は Tellus OS のタイル座標を表示したときの `Z` に相当します

![スクリーンショット 2022-11-15 16.48.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b47601f1-3b4b-ad04-c4a3-836c5fd89368.png)

PNG 画像を取得する際、タイル座標の `Z` 、 `X` 、 `Y` を指定する必要があります

シーン情報からは緯度経度しか取得できないため、緯度経度とタイル座標の変換を行います

ただし、 Tellus OS の地図はメルカトル図法になっているため、単純に計算できません

以下の関数を定義して計算します

※参考記事の Python 実装を Elixir に翻訳しました

```elixir
lat_lon_to_tile = fn zoom, lat, lon ->
  tile_count = 2 ** zoom

  x_tile =
    (lon + 180) / 360 * tile_count
    |> floor()

  y_tile =
    lat * :math.pi() / 180
    |> then(&(:math.tan(&1) + (1 / :math.cos(&1))))
    |> :math.log()
    |> then(&(1 - &1 / :math.pi()) / 2 * tile_count)
    |> floor()

  {x_tile, y_tile}
end
```

SAR データの中心を含むタイル座標を取得します

```elixir
lat = (scenes_info["data"]["left_bottom_lat"] + scenes_info["data"]["right_top_lat"])/ 2
lon = (scenes_info["data"]["left_bottom_lon"] + scenes_info["data"]["right_top_lon"]) / 2

{x, y} = lat_lon_to_tile.(z, lat, lon)
```

![スクリーンショット 2022-11-15 17.11.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f6a41d67-fb24-fb94-6e3e-d246d312ee4d.png)

PNG 画像には何種類かあるため、関数でダウンロード処理を定義しておきます

```elixir
get_png = fn type, work_id, z, x, y ->
  auth_header = get_auth_header.(sar_product_id)

  png_url = "#{base_url}/api/v2/works/#{work_id}/pngs/#{type}s/#{z}/#{x}/#{y}.png"

  %HTTPoison.Response{body: body} =
    png_url
    |> HTTPoison.get!([auth_header])

  img_path = "sar_#{type}_#{work_id}_#{z}_#{x}_#{y}.png"

  File.write!(img_path, body)

  img_path
  |> Evision.imread()
  |> Kino.render()

  img_path
end
```

まず、 `before` （差分干渉ペアのうち、日時が前のもの）を取得します

```elixir
before_img_path = get_png.("before", work_id, z, x, y)
```

![スクリーンショット 2022-11-15 17.12.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/6f89da10-8a9e-3b21-c9e7-a015bb3d752d.png)

`after` （差分干渉ペアのうち、日時が後のもの）を取得します

```elixir
after_img_path = get_png.("after", work_id, z, x, y)
```

![スクリーンショット 2022-11-15 17.13.23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b988a8a5-d540-9af5-fdd0-0986a4daa509.png)

コヒーレンス（before と after の重なり具合）を取得します

```elixir
coherence_img_path = get_png.("coherence", work_id, z, x, y)
```

![スクリーンショット 2022-11-15 17.15.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/25b5540c-4cb1-0088-43ad-bff3de4eb3ee.png)

差分干渉を取得します

```elixir
diff_img_path = get_png.("fringe_diff", work_id, z, x, y)
```

![スクリーンショット 2022-11-15 17.16.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/fc9f8bcc-8de0-122e-02ab-41d1b95ca42a.png)

## TIFF画像の取得

PNG ではタイルに区切られたデータになっているので、 TIFF で全体の差分干渉データを取得することもできます

```elixir
get_tif = fn type, work_id ->
  auth_header = get_auth_header.(sar_product_id)

  tiff_url = "#{base_url}/api/v2/works/#{work_id}/tifs/#{type}.tif"

  %HTTPoison.Response{body: body} =
    tiff_url
    |> HTTPoison.get!([auth_header])

  tiff_path = "sar_#{type}_#{work_id}.tif"

  File.write!(tiff_path, body)

  tiff_path
  |> Evision.imread()
  |> Kino.render()

  tiff_path
end
```

```elixir
get_tif.("fringe_diff", work_id)
```

![スクリーンショット 2022-11-15 17.17.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/309db682-e207-cdd5-b723-eb057f15d643.png)

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

SAR 画像を地図に重ねてみます

白い透過画像だと見えにくいので、白を赤に、透過を黒に変換します

```
get_red_image = fn img_path ->
  r =
    Evision.imread(img_path)
    |> Evision.Mat.to_nx(Nx.BinaryBackend)
    |> Nx.slice([0, 0, 0], [256, 256, 1])
    |> Nx.transpose(axes: [2, 0, 1])

  bg =
    0
    |> Nx.broadcast({256, 256, 2})
    |> Nx.transpose(axes: [2, 0, 1])

  [bg, r]
  |> Nx.concatenate()
  |> Nx.transpose(axes: [1, 2, 0])
  |> Nx.as_type({:f, 64})
  |> Evision.Mat.from_nx_2d()
end
```

```elixir
red_mat = get_red_image.(before_img_path)
```

![スクリーンショット 2022-11-15 17.19.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b7a8726f-55da-0d67-d056-c2023f0acee3.png)

MapLibre を使って画像を地図にオーバーレイする際、画像は URL で指定する必要があります

画像ファイルの保存先を LiveView から URL 指定で参照できなかったため、画像をデータURL（BASE64）にします

https://maplibre.org/

https://github.com/livebook-dev/maplibre

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

```elixir
r_img_base64 = get_data_url.(red_mat)
```

![スクリーンショット 2022-11-15 17.22.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/494b010e-b91a-d04d-eff0-f4d0c4a068b4.png)

MapLibre で画像を地図に重ねるとき、タイルの左右上下の緯度経度を指定する必要があります

なので、上で行っていたものと逆、タイル座標を緯度経度に変換する関数を定義します

```elixir
tile_to_lat_lon = fn zoom, x_tile, y_tile ->
  tile_count = 2 ** zoom

  lon_deg =
    x_tile / tile_count * 360 - 180

  lat_deg =
    :math.pi() * (1 - 2 * y_tile / tile_count)
    |> :math.sinh()
    |> :math.atan()
    |> Kernel.*(180 / :math.pi())

  {lat_deg, lon_deg}
end
```

```elixir
{bottom_lat, left_lon} = tile_to_lat_lon.(z, x, y)
```

![スクリーンショット 2022-11-15 17.30.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/28e16e84-0813-b7bb-2a6e-1f3eed99aac2.png)

```elixir
{top_lat, right_lon} = tile_to_lat_lon.(z, x + 1, y + 1)
```

![スクリーンショット 2022-11-15 17.48.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/470b96c6-803e-00e9-831c-5a33ea6af4e5.png)

MapLibre にエイリアスをつけておきます

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

地図表示用の関数を定義します

```elixir
show_map_overlay = fn z, {bottom_lat, left_lon, top_lat, right_lon}, img_base64 ->
  center = {(left_lon + right_lon) / 2, (bottom_lat + top_lat) / 2}

  # タイルの中心を地図の中心にする
  Ml.new(
    center: center,
    zoom: z,
    style: :terrain
  )
  # 画像をタイルの座標に配置する
  |> Ml.add_source(
    "sar-source",
    type: :image,
    url: img_base64,
    coordinates: [
      [left_lon, bottom_lat],
      [right_lon, bottom_lat],
      [right_lon, top_lat],
      [left_lon, top_lat],
    ]
  )
  # 画像をレイヤーとして地図に重ね、透過する
  |> Ml.add_layer(
    id: "overlay",
    source: "sar-source",
    type: :raster,
    layout: %{
      "visibility" => "visible"
    },
    paint: %{
      "raster-opacity" => 0.5
    }
  )
end
```

地図を表示します

```elixir
show_map_overlay.(10, {bottom_lat, left_lon, top_lat, right_lon}, r_img_base64)
```

![スクリーンショット 2022-11-15 17.49.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/fe49abe6-c37a-d130-2e02-d5733be0340a.png)

赤い SAR 画像を地図上に重ねることができました

他の衛星データも同じ要領で重ねることができますね

もっとズームしてみましょう

```elixir
z = 14
{x, y} = lat_lon_to_tile.(z, (bottom_lat + top_lat) / 2, (left_lon + right_lon) / 2)
```

![スクリーンショット 2022-11-15 17.36.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/15863316-6545-a1c1-ff67-144d42156cbd.png)

```elixir
before_img_path = get_png.("before", work_id, z, x, y)

r_img_base64 =
  before_img_path
  |> get_red_image.()
  |> get_data_url.()

{bottom_lat, left_lon} = tile_to_lat_lon.(z, x, y)
{top_lat, right_lon} = tile_to_lat_lon.(z, x + 1, y + 1)

show_map_overlay.(13, {bottom_lat, left_lon, top_lat, right_lon}, r_img_base64)
```

![スクリーンショット 2022-11-15 17.51.43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/38b0dfcf-63f7-b3c9-e43e-5c8035e667a5.png)

## まとめ

衛星データを画像として扱い、地図に重ねることができました

色々な可視化に応用できそうですね