Livebook で Tellus の標高データによる富士山とエベレストの画像・グラフ比較、海面上昇で沈む場所の可視化を実装する

## はじめに

日本発の衛星データプラットフォーム Tellus を使って標高データを取得します

https://www.tellusxdp.com/

今回は ASTER GDEM 3.0 という全世界の3次元地形データを使用します

ASTER GDEM はバージョン 3.0 から誰でも自由に利用可能になりました

Tellus からも無料でダウンロードできます

https://www.jspacesystems.or.jp/ersdac/GDEM/J/1.html

いつものように Elixir の Livebook を使っていきます

Livebook は簡単にデータを可視化することができます

https://livebook.dev/

実装したノートブックはこちら

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

## 出典

Credit: ASTER GDEM is courtesy of METI and NASA

## 事前準備

### 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 を呼ぶ準備ができました

## 実行環境

- Elixir: 1.14.2 OTP 24
- Livebook: 0.8.0

以下のリポジトリーの Docker コンテナ上で起動しています

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

:::note info

Docker が使える環境であれば簡単に実行できます

https://docs.docker.com/engine/install/

Docker Desktop を無償利用できない場合は Rancher Desktop を使ってください

https://rancherdesktop.io/

:::

私のリポジトリーを使う場合、以下の手順で Livebook を起動できます

```bash
git clone https://github.com/RyoWakabayashi/elixir-learning.git
```

```bash
cd elixir-learning
```

```bash
$ docker-compose up
...
Attaching to livebook
livebook    | [Livebook] Application running at http://0.0.0.0:<ポート番号>/?token=<認証トークン>
```

最後に表示される URL をブラウザで開けば Livebook にアクセスできます

![スクリーンショット 2022-09-30 17.40.59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8b51ac5a-f90a-2b99-f8e8-f8bac55f06ff.png)

右上 New notebook をクリックし、新しいノートブックを開けば準備完了です

![スクリーンショット 2022-09-30 17.43.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/85751035-145b-3e14-87af-7f7408069d14.png)

## セットアップ

必要なモジュールをインストールします

Livebook のセルに以下のコードを入力して実行してください

```elixir
Mix.install(
  [
    {:nx, "~> 0.4"},
    {:evision, "~> 0.1"},
    {:exla, "~> 0.4"},
    {:req, "~> 0.3"},
    {:kino, "~> 0.8"},
    {:kino_maplibre, "~> 0.1"},
    {:kino_vega_lite, "~> 0.1"}
  ],
  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クライアント
- [Kino](https://github.com/livebook-dev/kino): Livebook での入出力
- [Kino MapLibre](https://github.com/livebook-dev/kino_maplibre): 地図出力
- [Kino VegaLite](https://github.com/livebook-dev/kino_vega_lite): グラフ出力

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

## 情報の設定

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

```elixir
# Tellus のトークンを入力する
token_input = Kino.Input.password("Token")
```

## Tellus Traveler からデータを探す

データ取得には Tellus Satellite Data Traveler API を使用します

API 仕様は以下を参照

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

### データ取得用モジュール定義

データ取得用のモジュールを定義します

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

### データセットの確認

データセットの一覧を取得します

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

ダウンロード可能で説明文に「標高」が含まれるデータをテーブルに表示します

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

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

![データセット一覧](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/634cc972-2d9f-b7dc-ce92-3b40edca9711.png)

二つのデータセットが取得できました

それぞれ、説明文は以下のようになっています

- ASTER GDEM ver.3

> ASTER GDEM 3は、水平解像度30m相当、高さ方向精度約12m（日本域）の標高データです。
> 経済産業省が開発したASTERセンサのデータをもとに作成されており、値には建物や木の高さを含みます。
> 一部、データ欠損がある箇所はASTER以外のデータを用いており、建物や木の高さを含まない箇所もあります。
> パブリックドメインになりますので著作権表記は不要ですが、クレジットとして「ASTER GDEM is courtesy of METI and NASA」の表記を推奨します。

- AW3D30

> AW3D30は、水平解像度30m相当、高さ精度目標5mの標高データで高緯度地域もカバーされております。
> こちらのデータはJAXAの「だいち」(陸域観測技術衛星 ALOS)に搭載されたPRISMのデータをもとに作成されています。
> 
> ■ 解像度: 緯度・経度1秒 (30m相当)
> ■ 高さ精度：目標5m (標準偏差)
> ■ 関連URL：
> - プロダクト概要説明
> (https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d_j.htm)
> - オリジナルデータのフォーマット説明書(https://www.eorc.jaxa.jp/ALOS/jp/dataset/aw3d30/data/aw3d30v3.2_product_j_e1.2.pdf)※
> ※TellusOSや公開APIとはフォーマットが異なる場合があります。"

今回は ASTER GDEM の方を使うので、データセット ID に ASTER GDEM の ID を指定します

```elixir
dataset_id = "3f865d0b-6410-453f-b124-e0bf48544b45"
```

### 富士山のシーン選択

富士山のおおよその座標を指定します

```elixir
mt_fuji_location = {138.73, 35.36}
```

この座標をぐるっと囲む四角形を用意します

```elixir
mt_fuji_rectangle = [
  [
    [elem(mt_fuji_location, 0) - 0.001, elem(mt_fuji_location, 1) - 0.001],
    [elem(mt_fuji_location, 0) + 0.001, elem(mt_fuji_location, 1) - 0.001],
    [elem(mt_fuji_location, 0) + 0.001, elem(mt_fuji_location, 1) + 0.001],
    [elem(mt_fuji_location, 0) - 0.001, elem(mt_fuji_location, 1) + 0.001],
    [elem(mt_fuji_location, 0) - 0.001, elem(mt_fuji_location, 1) - 0.001]
  ]
]
```

Tellus から、この四角形を含む ASTER GDEM のシーン（ある場所、ある時間のデータ）の一覧を取得します

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

該当するのは1件だけです

GDEM の仕様を読むと分かりますが、 GDEM は世界中の陸地の座標を緯度経度1度ずつのタイルに分割して保有しています

日時も複数持ってはいないため、狭い範囲を指定すれば基本的には1件だけ該当します

https://www.jspacesystems.or.jp/ersdac/GDEM/J/2.html

```elixir
mt_fuji_scene_id =
  scenes_list
  |> Enum.at(0)
  |> then(& &1["id"])
```

これで `mt_fuji_scene_id` は `"61f5c40d-b0b6-42ec-8aca-dddeb964b534"` というシーン ID になりました

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

データセット ID とシーン ID を指定してデータをダウンロードします

```elixir
token_input
|> Kino.Input.read()
|> TellusTraveler.download(dataset_id, mt_fuji_scene_id)
```

実行結果は以下のように表示されます

```elixir
["/tmp/61f5c40d-b0b6-42ec-8aca-dddeb964b534/ASTGTMV003_N35E138_dem.tif",
 "/tmp/61f5c40d-b0b6-42ec-8aca-dddeb964b534/ASTGTMV003_N35E138_num.tif"]
```

GDEM の使用によると、 `..._dem.tif` は標高データ、 `..._num.tif` は QA （品質評価）データとなっています

`.tif` なので画像です

とりあえず Evision （OpenCV） で開いて表示してみましょう

一点気を付けないといけないのが、通常の画像と違い、色が符号付16bitになっている点です

GDEM は標高データなので、1m単位で表そうとすると、通常の色の範囲 0 〜 255 （符号なし8bit）では足りません

-65535 〜 65535 (符号付16bit) であれば海より低い場所もエベレストのように高い場所もカバーできます

`Evision.imread` する際、 `flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR()` を指定すれば 16bit 画像も問題なく開くことができます

ただし、表示上 255 より上は全て真っ白になります

```
"/tmp/#{mt_fuji_scene_id}"
|> File.ls!()
|> Enum.filter(fn filename -> Path.extname(filename) != ".txt" end)
|> Enum.sort()
|> Enum.map(fn filename ->
  ["/tmp", mt_fuji_scene_id, filename]
  |> Path.join()
  # 色が 16bit で格納されているため、 IMREAD_ANYDEPTH と IMREAD_ANYCOLOR を指定する
  |> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
  # 大きすぎるのでリサイズ
  |> Evision.resize({640, 640})
end)
|> Kino.Layout.grid(columns: 2)
```

![gdem.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/bf8a88ae-0868-ebcf-a269-4dadec1c7021.png)

左下が _dem.tif です

なんとなく右下が低いことは分かりますが、他は真っ白です

このままでは可視化できたとは言えないので分かりやすいように加工しましょう

まず _dem.tif を読み、テンソルにしておきます

```elixir
mt_fuji_dem =
  "/tmp/#{mt_fuji_scene_id}"
  |> File.ls!()
  |> Enum.find(fn filename -> String.ends_with?(filename, "_dem.tif") end)
  |> then(&Path.join(["/tmp", mt_fuji_scene_id, &1]))
  |> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
  |> Evision.Mat.to_nx(EXLA.Backend)
```

標高データの最大値と最小値を取り出します

```elixir
{mt_fuji_min_dig, mt_fuji_max_dig} = {
  mt_fuji_dem |> Nx.reduce_min() |> Nx.to_number(),
  mt_fuji_dem |> Nx.reduce_max() |> Nx.to_number()
}
```

結果は以下のようになります

```elixir
{-5, 3757}
```

最小値は -5 なので、一部海よりも低い場所があるようです

最大値は　3757 で、これが山頂です

誤差があるため実際の富士山頂とは少し違う値になっています

-5 から 3757 の範囲を 0 から 255 の範囲に変換します

```elixir
mt_fuji_dem_u8 =
  mt_fuji_dem
  |> Nx.subtract(mt_fuji_min_dig)
  |> Nx.multiply(255 / (mt_fuji_max_dig - mt_fuji_min_dig))
  |> Nx.as_type(:u8)

Evision.resize(mt_fuji_dem_u8, {640, 640})
```

![fuji_dem_gray.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f62b3c60-d0d7-cd7b-cefe-0553085a27a1.png)

山の稜線が見えました

白い筋が峰、黒い筋が谷です

そして一層白く輝いて見えるのが山頂です

ついでにカラーマップを適用して色を付けてみましょう

`COLORMAP_JET` は値の低い順に青 -> 水色 -> 緑 -> 黄 -> 赤と色が変化します

```elixir
mt_fuji_dem_color = Evision.applyColorMap(mt_fuji_dem_u8, Evision.Constant.cv_COLORMAP_JET())

Evision.resize(mt_fuji_dem_color, {640, 640})
```

![fuji_dem_color.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0dec3654-b7c9-9839-fa02-623b9fdb8bf3.png)

山頂は真っ赤でまるで噴火しているようですね

## 標高をグラフ化する

続いてグラフにしてみましょう

VegaLite は残念ながら 3D に対応していないため、地球を南北に両断したように、断面で標高を見てみます

断面標高グラフ表示用関数を用意します

```elixir
display_elevation_at_x = fn dem, x_index, width, max_y ->
  plot_data =
    dem[[x_index]]
    |> Nx.to_flat_list()
    |> Enum.with_index()
    |> Enum.map(fn {elevation, index} ->
      %{
        elevation: elevation,
        index: index
      }
    end)

  x_scale = %{"domain" => [0, 3601]}
  y_scale = %{"domain" => [-500, max_y]}

  VegaLite.new(width: width)
  |> VegaLite.data_from_values(plot_data)
  |> VegaLite.mark(:area)
  |> VegaLite.encode_field(:x, "index", type: :quantitative, scale: x_scale)
  |> VegaLite.encode_field(:y, "elevation", type: :quantitative, scale: y_scale)
end
```

画像左端、方位で言えば西端で切ってみましょう

```elixir
display_elevation_at_x.(mt_fuji_dem, 0, 700, 4000)
```

![fuji_graph_zero.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/f0ef10ec-07c4-daae-0d9d-be69f487e8fc.png)

こんな感じでグラフは山の形そのものになります

では最も高いところを見てみます

まず最も高い（標高の値が最も大きい）座標を探します

```elixir
max_index = mt_fuji_dem |> Nx.argmax() |> Nx.to_number()

max_x_index = div(max_index, 3601)
max_y_index = max_index - max_x_index * 3601

{max_x_index, max_y_index}
```

結果は以下のようになります

```elixir
{2302, 2618}
```

この座標の値を確認してみます

```elixir
mt_fuji_dem[[max_x_index, max_y_index]]
```

結果は以下のようになります

```elixir
#Nx.Tensor<
  s16
  EXLA.Backend<host:0, 0.2220939173.319946764.245111>
  3757
>
```

3757 なので、確かに先ほどの最大値と同じです

ではこの x 座標を指定して富士山をぶった斬りましょう

```elixir
display_elevation_at_x.(mt_fuji_dem, max_x_index, 700, 4000)
```

![fuji_graph_top.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/4e7f3dec-11b1-730f-4a1c-a0f1866ab95a.png)

なんと見事に尖っていることでしょう

もちろんグラフ全体の幅によって傾きは変わって見えますが、それにしてもよくこんな山に登る人が大勢いるものだと感心します

## グラフをアニメーションにする

一部だけを切り取ってみるのは勿体ないので、 3D ができない代わりにアニメーションで見てみましょう

@GeekMasahiro さんの記事を参考に、グラフを動かします

https://qiita.com/GeekMasahiro/items/5942bf1362f9caf8d16a

```elixir
x_scale = %{"domain" => [0, 3601]}
y_scale = %{"domain" => [-500, 4000]}

widget =
  VegaLite.new(width: 700)
  |> VegaLite.mark(:area)
  |> VegaLite.encode_field(:x, "index", type: :quantitative, scale: x_scale)
  |> VegaLite.encode_field(:y, "elevation", type: :quantitative, scale: y_scale)
  |> Kino.VegaLite.new()
```

まずウィジェットを用意すると、空のグラフが表示されます

続いて表示用関数を定義します

```elixiir
animate = fn dem, x_index ->
  plot_data =
    dem[[x_index]]
    |> Nx.to_flat_list()
    |> Enum.with_index()
    |> Enum.map(fn {elevation, y_index} ->
      %{
        elevation: elevation,
        index: y_index
      }
    end)

  Kino.VegaLite.clear(widget)
  Kino.VegaLite.push_many(widget, plot_data)
end
```

では 0 から 3600 まで値を動かしてみましょう

1 ずつだと凄く時間がかかるので、 20 ずつ動かします

```elixir
0..3600//20
|> Enum.map(fn x_index ->
  animate.(mt_fuji_dem, x_index)
  Process.sleep(100)
end)
```

![fuji_anime.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/a523b16a-95e5-a787-5ee3-b241eecf0da1.gif)

長いので記事には全編載せませんが（[また画像アップロードできなくなるので](https://qiita.com/RyoWakabayashi/items/3de067e1b1b5243a3762)）、山頂付近では山が隆起したように見えます

## 富士山と大分市を比較する

では富士山と我が町大分県大分市の標高を比較してみましょう

シーン ID 取得を関数化しておきます

```elixir
get_scene_id = fn location ->
  rectangle = [
    [
      [elem(location, 0) - 0.001, elem(location, 1) - 0.001],
      [elem(location, 0) + 0.001, elem(location, 1) - 0.001],
      [elem(location, 0) + 0.001, elem(location, 1) + 0.001],
      [elem(location, 0) - 0.001, elem(location, 1) + 0.001],
      [elem(location, 0) - 0.001, elem(location, 1) - 0.001]
    ]
  ]

  scenes_list =
    token_input
    |> Kino.Input.read()
    |> TellusTraveler.search(dataset_id, rectangle)

  scenes_list
  |> Enum.at(0)
  |> then(& &1["id"])
end
```

大分市のざっくりした座標を指定します

```elixir
oita_location = {131.64, 33.20}
```

大分市のシーン ID を取得します

```elixir
oita_scene_id = get_scene_id.(oita_location)
```

大分市の GDEM データをダウンロードします

```elixir
token_input
|> Kino.Input.read()
|> TellusTraveler.download(dataset_id, oita_scene_id)
```

比較するため、大分市も富士山と同じスケールでカラーマップに変換します

並べて表示してみましょう

```elixir
oita_dem =
  "/tmp/#{oita_scene_id}"
  |> File.ls!()
  |> Enum.find(fn filename -> String.ends_with?(filename, "_dem.tif") end)
  |> then(&Path.join(["/tmp", oita_scene_id, &1]))
  |> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
  |> Evision.Mat.to_nx(EXLA.Backend)

oita_dem_u8 =
  oita_dem
  |> Nx.subtract(mt_fuji_min_dig)
  |> Nx.multiply(255 / (mt_fuji_max_dig - mt_fuji_min_dig))
  |> Nx.as_type(:u8)

oita_dem_color = Evision.applyColorMap(oita_dem_u8, Evision.Constant.cv_COLORMAP_JET())

[mt_fuji_dem_color, oita_dem_color]
|> Kino.Layout.grid(columns: 2)
```

![fuji_oita.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/cc4ce0be-a6a2-553a-2c05-01d34ee34ce0.png)

左下の少し明るいところが九重ですね（これでもかなり高いはずですが）

https://travel.yahoo.co.jp/kanko/prefecture-oita/area-kokonoe/

大分も山だらけではありますが、いかに富士山が突き抜けて高いかが分かります

ではグラフを並べてみましょう

どちらも最も高いところです

```elixir
oita_max_index = oita_dem |> Nx.argmax() |> Nx.to_number()
oita_max_x_index = div(oita_max_index, 3601)

[
  display_elevation_at_x.(mt_fuji_dem, max_x_index, 300, 4000),
  display_elevation_at_x.(oita_dem, oita_max_x_index, 300, 4000)
]
|> Kino.Layout.grid(columns: 2)
```

![fuji_oita_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/1c2ca126-3d61-543d-cc1d-86dc440d9e29.png)

やはり富士山の尖り具合が強すぎる

こうして見ると地形の違いが直感的に分かりますね

## 富士山とエベレストを比較する

GDEM には世界中の陸地が入っているので、せっかくなら世界と比べてみましょう

というわけで文字通り世界最高峰のエベレストと比較してみます

```elixir
mt_everest_location = {86.92, 27.99}
```

```elixir
mt_everest_scene_id = get_scene_id.(mt_everest_location)
```

```elixir
token_input
|> Kino.Input.read()
|> TellusTraveler.download(dataset_id, mt_everest_scene_id)
```

富士山よりもエベレストの方がもちろん高いので、今度はエベレストのスケールに他を合わせます

```elixir
mt_everest_dem =
  "/tmp/#{mt_everest_scene_id}"
  |> File.ls!()
  |> Enum.find(fn filename -> String.ends_with?(filename, "_dem.tif") end)
  |> then(&Path.join(["/tmp", mt_everest_scene_id, &1]))
  |> Evision.imread(flags: Evision.Constant.cv_IMREAD_ANYDEPTH() + Evision.Constant.cv_IMREAD_ANYCOLOR())
  |> Evision.Mat.to_nx(EXLA.Backend)

{mt_everest_min_dig, mt_everest_max_dig} = {
  mt_everest_dem |> Nx.reduce_min() |> Nx.to_number(),
  mt_everest_dem |> Nx.reduce_max() |> Nx.to_number()
}
```

```elixir
min_dig = Enum.min([mt_fuji_min_dig, mt_everest_min_dig])
max_dig = Enum.max([mt_fuji_max_dig, mt_everest_max_dig])

{min_dig, max_dig}
```

ヒートマップ化するのは関数にしておきます

```elixir
get_heatmap = fn dem, min_dig, max_dig ->
  dem
  |> Nx.subtract(min_dig)
  |> Nx.multiply(256 / (max_dig - min_dig))
  |> Nx.as_type(:u8)
  |> Evision.applyColorMap(Evision.cv_COLORMAP_JET())
end
```

エベレスト、富士山、大分市を並べます

```elixir
[
  get_heatmap.(mt_everest_dem, min_dig, max_dig),
  get_heatmap.(mt_fuji_dem, min_dig, max_dig),
  get_heatmap.(oita_dem, min_dig, max_dig)
]
|> Kino.Layout.grid(columns: 3)
```

![everest.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/299653be-f1dd-6789-8b5e-5f75b4ef3086.png)

エベレストの山頂は右上ですね

山の高さも凄いですが、谷との落差が凄まじいです

エベレスト基準にすると、富士山もかなり大人しくなっています

ましてや大分は全体的に真っ暗、まるで海の中です

次にグラフを並べます

```elixir
mt_everest_max_x_index =
  mt_everest_dem
  |> Nx.argmax()
  |> Nx.to_number()
  |> div(3601)

[
  display_elevation_at_x.(mt_everest_dem, mt_everest_max_x_index, 200, 9000),
  display_elevation_at_x.(mt_fuji_dem, max_x_index, 200, 9000),
  display_elevation_at_x.(oita_dem, oita_max_x_index, 200, 9000)
]
|> Kino.Layout.grid(columns: 3)
```

![everest_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/87b26efc-69fd-a24a-694a-a5bcef896374.png)

エベレストはもう全体的に高すぎます

富士山が丸っと隠れてしまうような圧倒的高さです

しかも富士山以上に先端が尖っています

触れれば確実に刺さります

エベレストと比較してしまうと、大分県内は全域平野のようなものです

## 海面上昇シミュレーション

標高が分かるということは、例えば海面が現在より10m上がった場合、どこが沈むのかを見ることもできます

@zacky1972 さんが記事にしていますが、 2100 年には 10m 上がっていてもおかしくないのです

https://qiita.com/zacky1972/items/f4ed1bb4dd45f7d1e326#faq92-%E6%B5%B7%E9%9D%A2%E6%B0%B4%E4%BD%8D%E3%81%AF%E4%BB%8A%E5%BE%8C%E6%95%B0%E5%8D%81%E5%B9%B4%E9%96%93%E3%81%A7%E3%81%A9%E3%81%AE%E7%A8%8B%E5%BA%A6%E4%B8%8A%E6%98%87%E3%81%99%E3%82%8B%E3%81%8B

大分市で海に沈む場所を可視化してみましょう

標高が 0 （現在の海面) より大きく、 10 以下の場所を抽出します

```elixir
alpha =
  Nx.logical_and(Nx.greater(oita_dem, 0), Nx.less_equal(oita_dem, 10))
  |> Nx.select(255, 0)
  |> Nx.new_axis(-1)

alpha
|> Nx.as_type(:u8)
|> Evision.Mat.from_nx_2d()
|> Evision.resize({640, 640})
```

![oita_water.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/9e325a32-ebc9-7ea0-c085-8648882c2c58.png)

白いところが沈む場所です

海岸線に沿っているので、なんとなく大分県の形が見えますね

上の方に少し見えているのは山口県です

白だと見にくいのと、危機感を出すために赤くして他を透明にします

```elixir
bgra =
  [0, 0, 255]
  |> Nx.tensor()
  |> Nx.tile([3601, 3601, 1])
  |> then(&Nx.concatenate([&1, alpha], axis: 2))
  |> Nx.as_type(:u8)
  |> Evision.Mat.from_nx_2d()

Evision.resize(bgra, {640, 640})
```

![oita_red.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/078da249-34c0-def8-f81a-d6671a4a865b.png)

地図上に表示するため、画像をデータURLに変換する関数を用意します

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

MapLibre で地図上にこの画像を投影します

:::note warn
1 / 3600 上下に足し引きしているのは GDEM の仕様に沿っていますが、緯度方向に少しずれているのは恣意的に調整しています
:::

```elixir
center = {131.5, 33.5}

coordinates =
  [
    [131 - 1 / 3600 / 2, 34 - 9 / 3600 / 2],
    [132 + 1 / 3600 / 2, 34 - 9 / 3600 / 2],
    [132 + 1 / 3600 / 2, 33 - 11 / 3600 / 2],
    [131 - 1 / 3600 / 2, 33 - 11 / 3600 / 2],
  ]

img_base64 = get_data_url.(bgra)

MapLibre.new(center: center, zoom: 8, style: :terrain)
|> MapLibre.add_source("image", type: :image, url: img_base64, coordinates: coordinates)
|> MapLibre.add_layer(id: "overlay", source: "image", type: :raster, paint: %{"raster-opacity" => 0.5})
```

![oita.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/ea76e762-3096-4102-e8de-8673a6c76662.png)

大分市は海に面しているのもあって、10mも海面が上昇するとかなり沈んでしまいます

そして、残念なことに弊社 oec も海の底です

![oec.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/eb58c841-e812-3ab3-0806-562be3ac664b.gif)

## まとめ

標高データを使うことで、地形を直感的に捉えることができました

また、海面上昇の危機を改めて実感しました

持続可能な社会のために、私たちができることをやっていきましょう

https://www.mofa.go.jp/mofaj/gaiko/oda/sdgs/index.html

私はとにかく Elixir を推していきます！

@zacky1972 さんの連載記事で、 Elixir で SDGs に貢献できることを示してくれています

https://qiita.com/zacky1972/items/536e93af8404ed63b382