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

はじめに

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

前回の記事

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

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

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

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

Original data provided by NEC

参考記事

実行環境

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

Tellus の使い方

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

TelluSAR　の購入（無料）

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

トークンの作成

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

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

セットアップ

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

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

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

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

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

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

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

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

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

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

satellite_id = "asnaro2"

Tellus の認証

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

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

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

シーン一覧の取得

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

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

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

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

シーン情報の取得

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

scene_id = "AS200534028728-190103"

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

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

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

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

差分干渉データ作成依頼

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

after_scene_id = "AS200555328728-190117"

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

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

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

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

work_id = work["data"]["work_id"]

処理状況を確認します

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 になっています

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

PNGの取得

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

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

z = 11

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

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

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

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

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

※参考記事の Python 実装を 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 データの中心を含むタイル座標を取得します

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)

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

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 （差分干渉ペアのうち、日時が前のもの）を取得します

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

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

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

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

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

差分干渉を取得します

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

TIFF画像の取得

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

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

get_tif.("fringe_diff", work_id)

地図へのオーバーレイ

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

red_mat = get_red_image.(before_img_path)

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

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

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

r_img_base64 = get_data_url.(red_mat)

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

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

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

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

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

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

alias MapLibre, as: Ml

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

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

地図を表示します

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

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

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

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

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

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)

まとめ

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

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