衛星データプラットフォーム Tellus と Elixir Livebook を使って、植物の分布から田畑の利用状況を推測する（データ取得編）

## はじめに

経済産業省の「衛星データ利用環境整備・ソリューション開発支援事業」において、弊社は衛星データ無料利用事業者として支援を受けています

https://sdu.go.jp/

本事業では日本発の衛星データプラットフォーム「Tellus」のデータの一部を無料利用できます

（一般利用者が無料利用できるデータもあります）

https://www.tellusxdp.com/

「Tellus」から取得できる様々なデータを使って何をどう実装できるのか、現在調査検証中です

連載で「Tellus」の公式メディア「宙畑-sorabatake-」の記事を参考にしながら、 Elixir で衛星データを扱っています

https://sorabatake.jp/

参考にする記事

https://sorabatake.jp/28126/

今回も Elixir の Livebook を使っていくので、まだ読んでいない方は連載の環境構築編をご一読ください

連載記事

- [環境構築編](https://qiita.com/RyoWakabayashi/items/bf0590b5686c7462b4a0)
- データ取得編 <- ここ
- [NDVI算出編](https://qiita.com/RyoWakabayashi/items/b4fa2f7efa2a5d249a44)
- [市区町村切り抜き編](https://qiita.com/RyoWakabayashi/items/d65dd3a508284e67c8cd)
- [休耕地探索編](https://qiita.com/RyoWakabayashi/items/4039649daa0d22a3bda0)

本記事で実装したノートブックはこちら

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

## 出典

### 行政区域データ

「国土数値情報（行政区域データ）」（国土交通省）（<https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html>）を加工して作成

### 衛星データ

提供：だいち(ALOS)　AVNIR-2　データ(JAXA)

## Tellus のアカウント作成

Tellus から衛星データを取得するため、まずは Tellus のアカウントを作ります

公式ガイドを参照してください

https://www.tellusxdp.com/ja/market-service/service/guide/user-profile/20200220_000189.html

## トークンの作成

Tellus の API を呼び出すためのトークンを作成します

画面右上のログインユーザー名->「開発環境」をクリックします

![スクリーンショット 2022-11-11 18.03.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/b9d8cc19-6ccb-b712-033b-8715de3974ec.png)

「トークンの発行」をクリックすると、新しいトークンが作成され、下に表示されます

![スクリーンショット 2022-11-11 18.03.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/1ec22410-9b11-b73b-d43e-16c7cdbccdcf.png)

これで API を呼ぶ準備ができました

## 衛星データの選び方（Tellus Traveler）

[参考元の記事](https://sorabatake.jp/28126/)では、とりあえず AVNIR-2 のデータを使うことが前提になっていますが、そもそもどうやってデータを選べばいいのでしょうか

Tellus には Tellus Traveler という機能があり、ここでどのようなデータがあるのか探索することができます

https://www.tellusxdp.com/traveler/viewer/

ここでデータセットを選択したり、データの観測日や地理情報を指定することで検索が可能です

「どのデータセットか？」というのがよく分からなくても、おそらく「どこのデータか？」は決まっていると思います

例えば参考元の記事で使っている青森県藤崎町近辺を指定して検索してみましょう

![スクリーンショット 2023-01-05 18.16.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7220c4da-2613-89c7-6f7a-1ab0a6fe36e7.png)

すると、検索結果はこうなります

![スクリーンショット 2023-01-05 18.15.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/d59ac8cf-a667-c65d-2131-63c241e41f4f.png)

青い四角が該当する衛星データの領域なのですが、重なりすぎていて何が何だか分かりません

場所だけを指定してもデータセットを選択していないと、該当件数が多すぎて全く絞れていないのです

ではデータセットを「【Tellus公式】AVNIR-2_1B1」だけに絞ってみましょう

![スクリーンショット 2023-01-05 18.24.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/2f95e11d-f18e-7cef-3cfe-6393bf86572c.png)

そうするとかなり絞り込まれますが、これでも枠が重なり過ぎていてよく分かりません

結局一件一件詳細を開かないと、使えるデータなのかどうか確かめるのは難しいです

しかもここで恣意的に選んだデータの ID をプログラムに直書きするような形になるのはあまり美しいとは言えません

というわけで、 Elixir Livebook でこの問題を解決し、明確に衛星データを選択してみます

## Elixir Livebook による衛星データの選別

### セットアップ

Livebook を起動し、新しいノートブック（もしくは私の実装したノートブック）を開いてください

一番先頭のセル（セットアップセル:「Notebook dependencies and setup」 と書いてある場所）に以下の内容を入力してください

キーボードの「Shift+Enter」でセルの内容を実行します

```elixir
Mix.install(
  [
    {:nx, "~> 0.4"},
    {:evision, "~> 0.1"},
    {:exla, "~> 0.4"},
    {:req, "~> 0.3"},
    {:geo, "~> 3.4"},
    {:kino, "~> 0.8"},
    {:kino_maplibre, "~> 0.1.3"}
  ],
  config: [
    nx: [
      default_backend: EXLA.Backend
    ]
  ]
)
```

ここでは以下の依存モジュールをインストールしています

- [Nx](https://github.com/elixir-nx/nx/tree/main/nx): 行列演算
- [Evision](https://github.com/cocoa-xu/evision): 画像処理
- [EXLA](https://github.com/elixir-nx/nx/tree/main/exla): 行列演算の高速化
- [Req](https://github.com/wojtekmach/req): HTTPクライアント
- [Geo](https://github.com/bryanjos/geo): 地理情報システム
- [Kino](https://github.com/livebook-dev/kino): Livebook での入出力
- [Kino MapLibre](https://github.com/livebook-dev/kino_maplibre): 地図出力

また、 Nx による行列演算を EXLA を使って高速化する設定をしています

### 認証情報の設定

実行後、下のセル（黒いブロック）に以下の内容を入力し、実行します

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

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

### モジュールの作成

次のセルに以下の内容を入力します

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

Tellus Traveler の API を呼び出すためのモジュールです

API 仕様はこちら

https://www.tellusxdp.com/docs/travelers/

単に API へ HTTP リクエストを投げているだけですが、ひとまず内容は気にせず実行しましょう

次に作成したモジュールを使って、データセットの一覧を取得します

お金を使いたくないので、無償利用可能なデータを取得してみます

```elixir
datasets =
  token_input
  |> Kino.Input.read()
  |> TellusTraveler.get_datasets(false)
```

実行すると、以下のような結果が表示されます

![スクリーンショット 2023-01-05 18.41.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/819d9185-1a64-f4a9-3a44-835dc78610ee.png)

Tellus Traveler の API から返ってきた JSON 形式のデータです

このままでは見にくいので、必要な項目に絞ってテーブル表示します

```elixir
datasets
|> Enum.map(fn dataset ->
  %{
    "name" => dataset["name"],
    "description" => dataset["description"],
    "allowed_in" => dataset["permission"]["allow_network_type"]
  }
end)
|> Kino.DataTable.new(keys: ["name", "description", "allowed_in"])
```

![スクリーンショット 2023-01-05 18.43.26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/572a3b25-afb7-15d8-646c-d43da738bf94.png)

34 種類、無償で使えるデータセットがありました

ちなみにこのテーブルは並び替えや検索が可能で、そのまま Excel などに貼り付けることもできます

https://qiita.com/RyoWakabayashi/items/faf5c9dfd4a75d35daa6

### データセットの選び方

`allowed_in` の項目は「どこで使えるか」を表します

- tellus: Tellus 内のみでダウンロード不可
- global: ダウンロード可能

データを加工したり分析したりしたいので、 global のデータセットを選択しましょう

また、 `description` の項目にデータセットの概要が書かれています

今回は植生=植物がどれくらい生えているかを知りたいので、その場合には赤色光および近赤外線のデータが必要です

これらのデータは「光学」データなので、 `description` に「光学」が含まれるデータセットに絞ります

```elixir
opt_datasets =
  datasets
  |> Enum.filter(fn dataset ->
    dataset["permission"]["allow_network_type"] == "global" &&
    String.contains?(dataset["description"], "光学")
  end)

Kino.DataTable.new(opt_datasets, keys: ["name", "description"])
```

![スクリーンショット 2023-01-05 18.49.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f87fa70e-6f57-84c0-cd6b-3906e2b547b0.png)

これで4種類まで絞り込まれました

ちなみに、光学ではなく SAR 衛星のデータが欲しい場合は以下のようなコードになります

```elixir
sar_datasets =
  datasets
  |> Enum.filter(fn dataset ->
    dataset["permission"]["allow_network_type"] == "global" &&
    String.contains?(dataset["description"], "SAR")
  end)

Kino.DataTable.new(sar_datasets, keys: ["name", "description"])
```

### シーンの選び方

#### シーン一覧取得

あるデータセットの、ある日時の、ある領域の衛星データを「シーン」と読んでいます

今回は青森県藤崎町の領域を含む「シーン」を選択していきます

まず、データセットは先ほど絞り込んだ無償の光学データ4種類を対象にします

```elixir
dataset_id_list = Enum.map(opt_datasets, & &1["id"])
```

続いて、青森県藤崎町の地理情報（緯度経度で地球上のどの範囲にあるのか）を指定したいので、国土交通省の行政区域データを取得してきます

取得方法や詳細な利用方法は以下の記事を参照してください

https://qiita.com/RyoWakabayashi/items/1a12439d2971444ba33b

今回は青森県の行政区域データをダウンロードし、 /tmp/GML 内に展開しておきます

データを配置したディレクトリーを指定します

```elixir
gml_dir = "/tmp/GML"
```

データのうち、 GeoJSON ファイルを読み込みます

```elixir
geojson_file =
  gml_dir
  # ファイル一覧取得
  |> File.ls!()
  # `.geojson` で終わるもののうち先頭を取得
  |> Enum.find(&String.ends_with?(&1, ".geojson"))
```

```elixir
geojson_data =
  [gml_dir, geojson_file]
  |> Path.join()
  |> File.read!()
  |> Jason.decode!()
  |> Geo.JSON.decode!()
```

藤崎町（行政コード = 02361）だけを抜き出します

行政コードは以下のサイトを参照してください

https://nlftp.mlit.go.jp/ksj/gml/codelist/AdminAreaCd.html

```elixir
city_geojson = %Geo.GeometryCollection{
  geometries: Enum.filter(geojson_data.geometries, &(&1.properties["N03_007"] == "02361"))
}
```

とりあえず藤崎町を地図上にプロットしてみましょう

```elixir
city_map =
  MapLibre.new(center: {140.5, 40.7}, zoom: 7)
  |> MapLibre.add_geo_source("city_geojson", city_geojson)
  |> MapLibre.add_layer(
    id: "city_geojson",
    source: "city_geojson",
    type: :fill,
    paint: [fill_color: "#00ff00", fill_opacity: 0.5]
  )
```

![スクリーンショット 2023-01-05 19.07.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f0b17a4f-07a4-2b9e-0c96-8ec70b74e4fa.png)

ちゃんと青森県内に藤崎町らしき形が緑色でプロットされています

Tellus Traveler API では、多角形（Polygon）で領域を指定することで、その領域を含むシーンを検索することができます

ただし、指定できる多角形の座標数は 20 までになっているため、行政区域データのようなかなり角の多い多角形は指定できません

行政区域データの緯度・経度の最大値・最小値をそれぞれ計算し、藤崎町に外接する四角形を計算しましょう

```elixir
city_polygons =
  city_geojson.geometries
  |> Enum.at(0)
  |> then(& &1.coordinates)
```

```elixir
longitudes =
  city_polygons
  |> List.flatten()
  |> Enum.map(&elem(&1, 0))

latitudes =
  city_polygons
  |> List.flatten()
  |> Enum.map(&elem(&1, 1))

bbox = %{
  min_longitude: Enum.min(longitudes),
  max_longitude: Enum.max(longitudes),
  min_latitude: Enum.min(latitudes),
  max_latitude: Enum.max(latitudes),
}
```

```elixir
city_rectangle =
  [
    [
      [bbox.min_longitude, bbox.min_latitude],
      [bbox.max_longitude, bbox.min_latitude],
      [bbox.max_longitude, bbox.max_latitude],
      [bbox.min_longitude, bbox.max_latitude],
      [bbox.min_longitude, bbox.min_latitude],
    ]
  ]
```

四角形ですが、 Polygon では一周して元の位置に戻らないといけないので、最初と最後は同じ点で合計5ヶ所の緯度経度を指定しています

この四角形を使ってシーンを検索します

```elixir
scenes_list =
  token_input
  |> Kino.Input.read()
  |> TellusTraveler.search(dataset_id_list, city_rectangle)
```

![スクリーンショット 2023-01-05 19.16.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7dd99967-2d2e-0f04-d520-9fcf24d81349.png)

シーン一覧が取得できました

#### シーン領域の可視化

シーン毎に衛星データの領域を四角形で持っています

この領域を地図上にプロットしてみましょう

```elixir
scenes_list
|> Enum.map(& &1["geometry"]["coordinates"])
|> Enum.map(fn coordinates ->
  coordinates
  |> Enum.at(0)
  |> Enum.map(& List.to_tuple(&1))
end)
|> then(fn coordinates ->
  %Geo.GeometryCollection{
    geometries: [%Geo.Polygon{
      coordinates: coordinates
    }]
  }
end)
|> then(fn geojson ->
  city_map
  |> MapLibre.add_geo_source("area", geojson)
  |> MapLibre.add_layer(
    id: "area",
    source: "area",
    type: :line,
    paint: [line_color: "#00ff00"]
  )
end)
```

![スクリーンショット 2023-01-05 20.44.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8aa27255-a185-563b-5d75-b630407a457b.png)

藤崎市の領域を含むシーンといっても、どうやらそれぞれ領域が異なっているようです

しかもおおよそ4種類の領域に分かれていそうです

領域によっては藤崎市にかすっているだけで、ほとんど使えなさそうなものもあるように見えます

まだまだ絞り込みが必要です

#### データセット毎の件数確認

データセット毎に何件取得できたか確認してみます

```elixir
target_dataset_id_list =
  scenes_list
  # データセットIDでグルーピング
  |> Enum.group_by(& &1["dataset_id"])
  # それぞれの件数を取得
  |> Enum.map(fn {key, value} -> {key, Enum.count(value)}end)
  |> Enum.into(%{})
```

![スクリーンショット 2023-01-05 19.18.59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/45d20e17-2b25-6b89-d5af-356aa0157774.png)

4種類データセットを指定していましたが、結局藤崎町の領域にデータがあるのは1種類だけでした

もしここで複数データセットがある場合は複数データセットを横断して使うのか、最もシーン数が多いデータセットを使うのかなど検討してください

データセットIDからデータセット名や概要を取得してみます

```elixir
target_dataset_id_list
|> Enum.map(fn {dataset_id, count} ->
  Enum.find(opt_datasets, & &1["id"] == dataset_id)
  |> Map.merge(%{"count" => count})
end)
|> Kino.DataTable.new(keys: ["id", "name", "count", "description"])
```

![スクリーンショット 2023-01-05 19.21.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/c5ecdafc-bf04-73bc-d493-bf93c84427cf.png)

というわけで選ばれたのは「【Tellus公式】AVNIR-2_1B1」でした

参考元の記事と同じものです

これ以降、このデータセットだけを対象として扱います

```elixir
target_dataset_id = "ea71ef6e-9569-49fc-be16-ba98d876fb73"
```

#### データセットと雲量によるシーンの絞り込み

データセットIDと雲量を使って絞り込みます

`eo:cloud_cover` は 0 から 100 の範囲で雲の量を表します

この値が大きければ空が雲で覆われていて地面が観測できていない、ということになります

なので `eo:cloud_cover` に 25 未満を指定します

```elixir
target_scenes_list =
  scenes_list
  |> Enum.filter(fn scene ->
    scene["dataset_id"] == target_dataset_id &&
    scene["properties"]["eo:cloud_cover"] < 25
  end)
```

この時点で 31 件まで絞られます

#### 領域によるシーンの絞り込み

続いて領域（緯度経度）で絞り込みたいと思います

同じデータセットであれば衛星で撮影できる面積は同じなので、左上座標だけを使ってグループ分けしてみましょう

Evision (OpenCV) を使って「k-平均法」によるグループ分けを行います

各シーンから左上座標を取得し、テンソルに変換します

```elixir
coordinates_tensor =
  target_scenes_list
  |> Enum.map(fn scene ->
    scene["geometry"]["coordinates"]
    |> Enum.at(0)
    |> Enum.at(0)
  end)
  |> Nx.tensor()
```

![スクリーンショット 2023-01-05 19.32.05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/6c2466c4-06b5-89a6-522f-97aa195dfbed.png)

これを `Evision.kmeans` でグループ分けします

先ほど全シーンをプロットしたとき、おおよそ 4 種類の領域がありそうだったので、グループ数を 4 にしています

その他のパラメータ詳細については OpenCV の説明を参照してください

http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_ml/py_kmeans/py_kmeans_opencv/py_kmeans_opencv.html

```elixir
labels =
  coordinates_tensor
  |> Evision.kmeans(
    # ４グループに分ける
    4,
    # 必要ないが nil 指定できないので適当なテンソルを指定する
    Nx.tensor([0], type: :f32),
    # 3 = TERM_CRITERIA_EPS(1) + TERM_CRITERIA_MAX_ITER(2)
    {3, 10, 1.0},
    # 試行回数
    10,
    # 中心初期化手法を使用
    Evision.Constant.cv_KMEANS_PP_CENTERS
  )
  |> then(fn {compactness, labels, _} ->
    IO.inspect("compactness: #{compactness}")

    labels
    |> Evision.Mat.to_nx()
    |> Nx.to_flat_list()
  end)
```

デバッグ出力している compactness が 0.01 くらいならおおよそグループ分けできています

compactness が大きよようであればグループ数を増やしてみましょう

もっと適切なグループ数を探したい場合はエルボー法を使いますが、それは後日別記事にします

![スクリーンショット 2023-01-05 19.48.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/63c92155-fa62-5231-e2fa-426057427bc7.png)

何はともあれ４グループに分けることができました

出力される List の 0 1 2 3 がグループの番号になっています

このグループ番号 List を使って、シーンをグループに分けます

```elixir
target_scenes_group =
  [target_scenes_list, labels]
  |> Enum.zip()
  |> Enum.group_by(fn {_, label} -> label end, fn {scene, _} -> scene end)
```

![スクリーンショット 2023-01-05 19.50.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/bcf62d55-49c3-9922-9cec-8f46d9234628.png)

グループ毎の件数を見てみます

```elixir
target_scenes_group
|> Enum.map(fn {label, scenes} -> {label, Enum.count(scenes)} end)
|> Enum.into(%{})
```

![スクリーンショット 2023-01-05 21.05.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/41a9f456-56f6-4a4c-dc38-920a55c9f7d7.png)

各グループ同じくらいの件数です

これだけではピンとこないので、各グループの先頭シーンを地図上にプロットしてみます

```elixir
# グループ毎に先頭シーンの領域を取得
scene_geojson_map =
  target_scenes_group
  |> Enum.map(fn {label, scene_list} ->
    scene_list
    |> Enum.at(0)
    |> then(& &1["geometry"]["coordinates"])
    |> Enum.at(0)
    |> Enum.map(& List.to_tuple(&1))
    |> then(fn coordinates ->
      %Geo.GeometryCollection{
        geometries: [%Geo.Polygon{
          coordinates: [coordinates]
        }]
      }
    end)
    |> then(fn geojson ->
      {label, geojson}
    end)
  end)
  |> Enum.into(%{})
```

```elixir
# 藤崎町の中心座標
center =
  {
    (bbox.min_longitude + bbox.max_longitude) / 2,
    (bbox.min_latitude + bbox.max_latitude) / 2
  }
```

まずシーンと藤崎町がどのように重なっているか見るため、地図上で藤崎町を黒く塗りつぶしておきます

```elixir
city_map =
  MapLibre.new(center: center, zoom: 7)
  |> MapLibre.add_geo_source("city_geojson", city_geojson)
  |> MapLibre.add_layer(
    id: "city",
    source: "city_geojson",
    type: :fill,
    paint: [fill_color: "#000000"]
  )
```

![スクリーンショット 2023-01-05 19.56.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/00567870-1075-cae2-35b3-5f09c8e2378a.png)

そして、ここからがこの記事のハイライトです

各グループの領域を地図にプロットし、タブで切り替えて表示できるようにします

```elixir
scene_geojson_map
|> Enum.map(fn {label, scene_geojson} ->
  map =
    city_map
    |> MapLibre.add_geo_source("area", scene_geojson)
    |> MapLibre.add_layer(
      id: "area",
      source: "area",
      type: :fill,
      paint: [fill_color: "#00ff00", fill_opacity: 0.5]
    )

  {Integer.to_string(label), map}
end)
|> Kino.Layout.tabs()
```

![tabs.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/48424066-55b2-2ae7-1782-aee6f01f5f81.gif)

グループ1は藤崎市がすみっこすぎて、一部がはみ出していそうです

グループ２は藤崎市をかすっているだけで使えそうにありません

グループ0が一番良さそうです

こんな可視化ができるのが Livebook のすごいところですね

Kino.Layout についてはこちら

https://qiita.com/RyoWakabayashi/items/b60450a15c3be021f1b5

https://qiita.com/RyoWakabayashi/items/fd69b0fc73d65b1e3df6

では、グループ内の各シーンが本当に藤崎町を含んでいるのかも確認してみましょう

```elixir
target_scenes_group[0]
|> Enum.map(fn scene ->
  coordinates =
    scene["geometry"]["coordinates"]
    |> Enum.at(0)
    |> Enum.map(& List.to_tuple(&1))

  scene_geojson =
    %Geo.GeometryCollection{
      geometries: [%Geo.Polygon{
        coordinates: [coordinates]
      }]
    }

  map =
    city_map
    |> MapLibre.add_geo_source("area", scene_geojson)
    |> MapLibre.add_layer(
      id: "area",
      source: "area",
      type: :fill,
      paint: [fill_color: "#00ff00", fill_opacity: 0.5]
    )

  # 観測日をタブにする
  date = String.slice(scene["properties"]["start_datetime"], 0..9)

  {date, map}
end)
|> Kino.Layout.tabs()
```

![tabs_2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/eb2cffb1-2aca-4872-b48b-3dd4a61859aa.gif)

どのシーンもほぼ同じ領域を示していますね

このグループのデータを使えば、同じ領域の推移を観測できそうです

シーンIDの一覧を見てみましょう

```elixir
scene_id_list =
  target_scenes_group[0]
  |> Enum.map(& &1["id"])
```

以下の結果が表示されます

```elixir
["202ce08d-ba4b-4ffe-8165-109fd3a8b917", "34d8dc6f-fdd1-4542-a038-c1235a5a97fa",
 "12ad308b-6ce1-40ec-9ebf-f0215c30500e", "e2e85b2e-a208-4a65-87fd-b92721b037a8",
 "208a3618-7700-421b-bf05-fd59551cc1aa", "d5ce7320-5b25-4ced-bda5-0e25a9d75940",
 "9d14706f-cee7-4eb4-9151-2558609c3de0", "3f4555ac-eaf3-4066-a1ba-20bb1ec1c0b3"]
```

このシーンIDは参考元記事で使ったノートブックで指定しているものとほぼ同じ（一つだけこちらが多い）です

https://github.com/sorabatake/softwaredesign_1/blob/master/Anarysis_Example_using_AVNIR2.ipynb

したがって、 Elixir Livebook を使うことで、 Tellus Traveler の画面検索と同等（もしくはより根拠を示せる形で）衛星データを選択できる、ということになります

## 衛星データのダウンロード

選択したシーンIDを使ってデータをダウンロードしましょう

データが大きいので、相当の時間が掛かるのと、ストレージが満杯にならないように気をつけてください

```elixir
scene_id_list
|> Enum.map(fn scene_id -> 
  token_input
  |> Kino.Input.read()
  |> TellusTraveler.download(target_dataset_id, scene_id)
end)
```

## ダウンロードした画像の表示

```elixir
scene_id = Enum.at(scene_id_list, 2)

"/tmp/#{scene_id}"
|> File.ls!()
|> Enum.filter(fn filename -> Path.extname(filename) != ".txt" end)
|> Enum.sort()
|> Enum.map(fn filename ->
  ["/tmp", scene_id, filename]
  |> Path.join()
  |> Evision.imread()
  # 大きすぎるのでリサイズ
  |> Evision.resize({640, 640})
end)
|> Kino.Layout.grid(columns: 2)
```

![スクリーンショット 2023-01-05 20.19.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/d19b87de-18e7-d11f-5cd3-2c9a06fffeb0.png)

十和田湖や陸奥湾が見えますね

下4枚の画像は各波長帯の光学データで、左下が赤色光、右下が赤外線のデータです

次回以降、これらの衛星データを使って植生を計算し、地図へのプロットや季節による推移のグラフ化を実装していきます

## まとめ

Elixir Livebook を使うことで、かなり視覚的、効果的に衛星データの選別を行うことができました

特にタブ表示はかなり強力です

今後、 Livebook がより広く使われることを期待します

次回は NDVI という値を計算することによって、植物の分布（植生）を可視化します

https://qiita.com/RyoWakabayashi/items/b4fa2f7efa2a5d249a44