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

## はじめに

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

https://sdu.go.jp/

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

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

https://www.tellusxdp.com/

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

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

https://sorabatake.jp/

参考にする記事

https://sorabatake.jp/28126/

今回は植生（植物の集まり、植物がどれくらい生えているか）を見るために算出した NDVI(Normalized Difference Vegetation Index) = 正規化植生指数について、特定の市区町村内だけを切り抜きます

これによって、市区町村毎の植生の比較や、市区町村内の植生推移が可視化できます

引き続き Elixir の Livebook を使っていきます

まだ読んでいない方は環境構築編、データ取得編から読んでください

連載記事

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

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

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

## 出典

### 行政区域データ

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

### 衛星データ

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

## 前提条件

[データ取得編](https://qiita.com/RyoWakabayashi/items/c87fe1256bcd105d370c) で対象としたシーンのデータを Tellus からダウンロードし、 `/tmp/<シーンID>/` 配下に保存しているものとします

## セットアップ

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

セットアップセルに以下の内容を入力して実行してください

```elixir
Mix.install(
  [
    {:nx, "~> 0.4"},
    {:evision, "~> 0.1"},
    {:exla, "~> 0.4"},
    {:geo, "~> 3.4"},
    {:jason, "~> 1.4"},
    {:kino, "~> 0.8"},
    {:kino_maplibre, "~> 0.1.3"},
    {:kino_vega_lite, "~> 0.1.7"}
  ],
  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): 行列演算の高速化
- [Geo](https://github.com/bryanjos/geo): 地理情報システム
- [Jason](https://github.com/michalmuskala/jason): JSONパーサー
- [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 を使って高速化する設定をしています

## データの指定

[データ取得編](https://qiita.com/RyoWakabayashi/items/c87fe1256bcd105d370c) で選択したシーンIDリストを指定しておきます

```elixir
scene_id_list = [
  "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"
]
```

## ヘッダーファイルの読込モジュール

[NDVI算出編](https://qiita.com/RyoWakabayashi/items/b4fa2f7efa2a5d249a44)で使ったヘッダーファイルの読込モジュールを用意します

```elixir
defmodule Coordinate do
  defstruct latitude: 0.0, longitude: 0.0
end

defmodule BandInfo do
  defstruct gain: 0.0, offset: 0.0
end

defmodule Affine do
  defstruct a: 0.0, b: 0.0, c: 0.0, d: 0.0
end

defmodule Conversion do
  defstruct x: 0.0, y: 0.0
end
```

```elixir
defmodule HeaderInfo do
  defstruct blue_band: %BandInfo{},
            green_band: %BandInfo{},
            red_band: %BandInfo{},
            nir_band: %BandInfo{},
            center: %Coordinate{},
            left_top: %Coordinate{},
            right_top: %Coordinate{},
            left_bottom: %Coordinate{},
            right_bottom: %Coordinate{},
            affine: %Affine{},
            conversion: %Conversion{},
            degree: 0.0,
            map_degree: 0.0,
            datetime: nil,
            product_id: ""

  def get_string(info, start, value) do
    info
    |> String.slice(start, value)
    |> String.trim()
  end

  def get_value(info, start, value) do
    info
    |> get_string(start, value)
    |> String.to_float()
  end

  def read(hdr_file_path) do
    info = File.read!(hdr_file_path)

    %HeaderInfo{
      # 青色光バンド
      blue_band: %BandInfo{
        gain: get_value(info, 1721, 8),
        offset: get_value(info, 1729, 8)
      },
      # 緑色光バンド
      green_band: %BandInfo{
        gain: get_value(info, 1737, 8),
        offset: get_value(info, 1745, 8)
      },
      # 赤色光バンド
      red_band: %BandInfo{
        gain: get_value(info, 1752, 8),
        offset: get_value(info, 1760, 8)
      },
      # 近赤外線バンド
      nir_band: %BandInfo{
        gain: get_value(info, 1768, 8),
        offset: get_value(info, 1776, 8)
      },
      # 画像中央
      center: %Coordinate{
        latitude: get_value(info, 248, 16),
        longitude: get_value(info, 264, 16)
      },
      # 画像左上
      left_top: %Coordinate{
        latitude: get_value(info, 376, 16),
        longitude: get_value(info, 392, 16)
      },
      # 画像右上
      right_top: %Coordinate{
        latitude: get_value(info, 408, 16),
        longitude: get_value(info, 424, 16)
      },
      # 画像左下
      left_bottom: %Coordinate{
        latitude: get_value(info, 440, 16),
        longitude: get_value(info, 456, 16)
      },
      # 画像右下
      right_bottom: %Coordinate{
        latitude: get_value(info, 472, 16),
        longitude: get_value(info, 488, 16)
      },
      affine: %Affine{
        a: get_value(info, 1224, 16),
        b: get_value(info, 1240, 16),
        c: get_value(info, 1256, 16),
        d: get_value(info, 1272, 16)
      },
      conversion: %Conversion{
        x: get_value(info, 1208, 8),
        y: get_value(info, 1216, 8)
      },
      degree: get_value(info, 760, 16),
      map_degree: get_value(info, 921, 16),
      datetime:
        info
        |> get_string(192, 24)
        |> then(fn str ->
          String.slice(str, 0, 4) <>
            "-" <>
            String.slice(str, 4, 2) <>
            "-" <>
            String.slice(str, 6, 2) <>
            "T" <>
            String.slice(str, 8, 2) <>
            ":" <>
            String.slice(str, 10, 2) <>
            ":" <> String.slice(str, 12, 2)
        end)
        |> NaiveDateTime.from_iso8601!(),
      product_id: get_string(info, 128, 16)
    }
  end
end
```

## NDVIの算出モジュール

[NDVI算出編](https://qiita.com/RyoWakabayashi/items/b4fa2f7efa2a5d249a44)で使ったNDVIの算出モジュールを用意します

```elixir
defmodule NDVI do
  def read_header(file_path_list) do
    file_path_list
    |> Enum.find(fn file -> Path.extname(file) == ".txt" end)
    |> HeaderInfo.read()
  end

  def get_band_tensor(file_path_list, prefix) do
    file_path_list
    |> Enum.find(fn file ->
      file
      |> Path.basename()
      |> String.starts_with?(prefix)
    end)
    |> Evision.imread(flags: Evision.Constant.cv_IMREAD_GRAYSCALE())
    |> Evision.Mat.to_nx(EXLA.Backend)
  end

  def calc(file_path_list) do
    header_info = read_header(file_path_list)

    red_tensor =
      file_path_list
      |> get_band_tensor("IMG-03")
      |> Nx.multiply(header_info.red_band.gain)
      |> Nx.add(header_info.red_band.offset)

    nir_tensor =
      file_path_list
      |> get_band_tensor("IMG-04")
      |> Nx.multiply(header_info.nir_band.gain)
      |> Nx.add(header_info.nir_band.offset)

    ndvi_tensor =
      Nx.select(
        Nx.multiply(
          Nx.not_equal(red_tensor, 0),
          Nx.not_equal(nir_tensor, 0)
        ),
        Nx.divide(
          Nx.subtract(nir_tensor, red_tensor),
          Nx.add(nir_tensor, red_tensor)
        ),
        0
      )

    {header_info, ndvi_tensor}
  end
end
```

各シーンの NDVI を算出します

```elixir
ndvi_list =
  scene_id_list
  |> Enum.map(fn scene_id ->
    "/tmp/#{scene_id}"
    |> File.ls!()
    |> Enum.map(fn filename -> Path.join(["/tmp", scene_id, filename]) end)
    |> NDVI.calc()
  end)
```

先頭のシーンについて NDVI を取得し、画像として表示します

```elixir
{header_info, ndvi_tensor} = Enum.at(ndvi_list, 0)

ndvi_map =
  ndvi_tensor
  |> Nx.multiply(128)
  |> Nx.add(128)
  |> Nx.as_type(:u8)
  |> Evision.Mat.from_nx_2d()
  |> Evision.applyColorMap(Evision.Constant.cv_COLORMAP_WINTER())

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

![スクリーンショット 2023-01-06 18.13.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/9169f4c0-df1c-bbc1-f83e-7650b87e1513.png)

これを地図に重ねてみます

```elixir
center = {header_info.center.longitude, header_info.center.latitude}
```

```elixir
coordinates = [
  [header_info.left_top.longitude, header_info.left_top.latitude],
  [header_info.right_top.longitude, header_info.right_top.latitude],
  [header_info.right_bottom.longitude, header_info.right_bottom.latitude],
  [header_info.left_bottom.longitude, header_info.left_bottom.latitude]
]
```

```elixir
# 地図にプロットするために BASE64 エンコードする
get_data_url = fn mat ->
  Evision.imencode(".png", mat)
  |> Base.encode64()
  |> then(&"data:image/png;base64,#{&1}")
end
```

```elixir
img_url =
  ndvi_map
  |> Evision.resize({4960, 4960})
  |> get_data_url.()

MapLibre.new(center: center, zoom: 8, style: :terrain)
|> MapLibre.add_source("image", type: :image, url: img_url, coordinates: coordinates)
|> MapLibre.add_layer(id: "overlay", source: "image", type: :raster)
```

![スクリーンショット 2023-01-06 18.15.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/7b466889-d2a0-8135-11d1-c0e6e5a4d72e.png)

ここまでは前回のおさらいです

## 市区町村のポリゴンデータを用意する

[国土交通省の行政区域データダウンロードページ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html)から青森県、令和4年のデータをダウンロードします

行政区域データについては以下の記事を参照してください

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

この記事ではダウンロードしたデータを `/tmp/GML` に展開したものとして進めます

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

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

行政コード = 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.68}, zoom: 11, style: :terrain)
  |> 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]
  )
```

![fujisaki.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/9d219075-90e9-f2fe-4f92-47ae27b8dac0.gif)

地図上の藤崎町の境界線と、緑色で表示した行政区域データが見事に重なっています

この行政区域データのポリゴン（多角形）で NDVI を切り抜き、藤崎町内の NDVI を可視化します

## NDVI 画像の補正

衛星データの形状を確認してみます

```elixir
{ndvi_height, ndvi_width, 3} = Evision.Mat.shape(ndvi_map)
```

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

```elixir
{8000, 7244, 3}
```

つまり、衛星データは縦 8,000 ピクセル × 横 7,244 ピクセル × 3 色の画像データになっています

しかし、先ほど地図上に NDVI をプロットしたのを見ても分かるとおり、この画像データは地軸に対して斜めに傾いています

どれくらい傾いているかはヘッダー情報に格納されています

```elixir
header_info.degree
```

値は以下のような度（°）で与えられます

```elixir
13.849887
```

この分回転して地軸に対して並行にすれば、行政区域データに重ねることができそうです

四隅の角が画像内に残るように画像を回転する関数を定義します

```elixir
wrap_rotate = fn mat, degree ->
  # 元の大きさ
  {src_height, src_width, 3} = mat.shape

  # 四隅がはみ出ることを考慮しない回転のアフィン変換
  affine =
    {src_width / 2, src_height / 2}
    |> Evision.getRotationMatrix2D(degree, 1)
    |> Evision.Mat.to_nx(EXLA.Backend)

  # 回転後に四隅を画像内に含むためのサイズを計算
  cos = Nx.abs(affine[[0, 0]])
  sin = Nx.abs(affine[[0, 1]])

  dst_width =
    Nx.add(Nx.multiply(src_height, sin), Nx.multiply(src_width, cos))
    |> Nx.to_number()
    |> trunc()

  dst_height =
    Nx.add(Nx.multiply(src_height, cos), Nx.multiply(src_width, sin))
    |> Nx.to_number()
    |> trunc()

  # (変換後サイズ - 変更前サイズ) / 2 の分だけずらして画像の中心を回転の中心に合わせる
  bias =
    Nx.tensor([
      [0, 0, (dst_width - src_width) / 2],
      [0, 0, (dst_height - src_height) / 2]
    ])

  # 四隅を画像内に含めるアフィン変換
  new_affine = Nx.add(affine, bias)

  Evision.warpAffine(mat, new_affine, {dst_width, dst_height})
end
```

NDVI の画像を回転させます

```elixir
rotated_map = wrap_rotate.(ndvi_map, -header_info.degree)

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

![rotated.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/25dba34c-5ffe-4fe8-dfc6-8b72f423ec1e.png)

また、画像は長方形になっていますが、緯度経度を見ると、これも少しずれていることが分かります

長方形であるならば、回転していたとしても上辺の経度の差と下辺の経度の差は同じはずです

```elixir
diff_longitude =
  (header_info.left_top.longitude - header_info.right_top.longitude) -
  (header_info.left_bottom.longitude - header_info.right_bottom.longitude)
```

しかし、この上辺と下辺の経度の差を計算すると、以下のような値になります

```elixir
-0.008156400000018493
```

この値の分だけ台形になっていそうです

正確には台形補正をすべきですが、ここでは単純に誤差の半分だけ左側に水平移動します

:::note warn
この辺りは正確ではないです

しかし何らかのか補正をしないとぴったり当てはまりません

衛星データに詳しい方、ご教示いただけると助かります
:::

1ピクセルあたりの緯度経度を計算します

```elixir
ndvi_longitudes = [
  header_info.left_top.longitude,
  header_info.right_top.longitude,
  header_info.left_bottom.longitude,
  header_info.right_bottom.longitude
]
```

```elixir
ndvi_latitudes = [
  header_info.left_top.latitude,
  header_info.right_top.latitude,
  header_info.left_bottom.latitude,
  header_info.right_bottom.latitude
]
```

```elixir
ndvi_min_longitude = Enum.min(ndvi_longitudes)
ndvi_max_longitude = Enum.max(ndvi_longitudes)
ndvi_min_latitude = Enum.min(ndvi_latitudes)
ndvi_max_latitude = Enum.max(ndvi_latitudes)

{ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude}
```

```elixir
{ndvi_height, ndvi_width, 3} = Evision.Mat.shape(rotated_map)
```

```elixir
pix_per_longitude = ndvi_width / (ndvi_max_longitude - ndvi_min_longitude)
pix_per_latitude = ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)

{pix_per_longitude, pix_per_latitude}
```

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

```elixir
{8483.693398865544, 11116.281262859115}
```

この値を使って、横方向への補正値を計算します

```elixir
diff_x = trunc(diff_longitude / 2 * pix_per_longitude)
```

補正値は以下の値になります

```elixir
-34
```

補正値を使って画像を横にずらします

```elixir
affine =
  [
    [1, 0, diff_x],
    [0, 1, 0]
  ]
  |> Nx.tensor()
  |> Nx.as_type(:f32)

corrected_map = Evision.warpAffine(rotated_map, affine, {ndvi_width, ndvi_height})

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

![corrected.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/0716b5fc-ab98-1c3c-2e61-9eb077ce31ef.png)

ほとんど変わっていませんが、僅かながら右上の角と右端の間に隙間ができています

その分、左下の角が画像からはみ出ました

## マスク画像による切り抜き

NDVI の画像を切り抜くためのマスク画像を作成します

```elixir
# 緯度経度の最大最小を求める
city_coordinates =
  city_geojson.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()

city_longitudes = Enum.map(city_coordinates, &elem(&1, 0))
city_latitudes = Enum.map(city_coordinates, &elem(&1, 1))

city_min_longitude = Enum.min(city_longitudes)
city_max_longitude = Enum.max(city_longitudes)
city_min_latitude = Enum.min(city_latitudes)
city_max_latitude = Enum.max(city_latitudes)

{city_min_longitude, city_max_longitude, city_min_latitude, city_max_latitude}
```

```elixir
# 最大最小の差からマスク画像の大きさを求める
city_width = trunc((city_max_longitude - city_min_longitude) * pix_per_longitude)
city_height = trunc((city_max_latitude - city_min_latitude) * pix_per_latitude)

{city_width, city_height}
```

行政区域データの各座標（緯度経度）をマスク画像内の座標（XY）に変換します

```elixir
# 緯度経度の最小から最大をピクセル数の0から幅・高さにスケールする
normalized_points =
  city_geojson.geometries
  |> Enum.map(& &1.coordinates)
  |> Enum.map(fn coordinate ->
    coordinate
    |> Enum.at(0)
    |> Enum.map(fn {x, y} ->
      [
        trunc((x - city_min_longitude) * pix_per_longitude),
        # 縦方向は緯度が北緯の場合逆転するため、高さから引く
        trunc(city_height - (y - city_min_latitude) * pix_per_latitude)
      ]
    end)
    |> Nx.tensor(type: :s32)
  end)
```

真っ黒で不透明な空画像を作ります

```elixir
# 空画像を用意する
empty_mat =
  [0, 0, 0, 255]
  |> Nx.tensor(type: :u8)
  |> Nx.tile([city_height, city_width, 1])
  |> Evision.Mat.from_nx_2d()

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

![black.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/409add73-13b5-fc85-c884-11e6016b8a62.png)

空画像内に行政区域データを多角形として書き込みます

```elixir
# ポリゴンを透明色で塗りつぶす
mask_mat = Evision.fillPoly(empty_mat, normalized_points, {0, 0, 0, 0})

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

![polygon.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/6c036bd0-2800-8544-09b6-d79b86dd7adc.png)

白い部分は透明になっています

## 地図へのプロット

マスク画像を地図にプロットしてみます

```elixir
city_center = {
  (city_min_longitude + city_max_longitude) / 2,
  (city_min_latitude + city_max_latitude) / 2
}
```

```elixir
city_bbox_coordinates =
  [
    [city_min_longitude, city_max_latitude],
    [city_max_longitude, city_max_latitude],
    [city_max_longitude, city_min_latitude],
    [city_min_longitude, city_min_latitude]
  ]
```

```elixir
image_url =
  mask_mat
  |> Evision.resize({trunc(city_height / 4), trunc(city_width / 4)})
  |> get_data_url.()

MapLibre.new(center: city_center, zoom: 11, style: :terrain)
|> MapLibre.add_source("city_mask", type: :image, url: image_url, coordinates: city_bbox_coordinates)
|> MapLibre.add_layer(id: "overlay", source: "city_mask", type: :raster)
```

![スクリーンショット 2023-01-11 17.15.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/75b52bc1-c338-3576-19db-d4289bf1af7a.png)

藤崎町をくり抜くことができています

## NDVI を行政区域で切り抜く

いよいよ NDVI を切り抜きます

まず、行政区域を含む領域を四角形で切り抜きます

```elixir
# 緯度経度を元に NDVI 画像内の行政区域の位置を計算する
city_left_top_x = trunc((city_min_longitude - ndvi_min_longitude) * pix_per_longitude)
city_left_top_y = trunc((ndvi_max_latitude - city_max_latitude) * pix_per_latitude)

{city_left_top_x, city_left_top_y}
```

```elixir
city_ndvi_bbox_map =
  corrected_map
  |> Evision.Mat.to_nx(EXLA.Backend)
  |> then(
    & &1[
      [
        city_left_top_y..(city_left_top_y + city_height - 1),
        city_left_top_x..(city_left_top_x + city_width - 1)
      ]
    ]
  )
  |> Evision.Mat.from_nx_2d()
```

![ndvi_city_bbox.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/069c5f25-697d-1d18-433c-d393b9bf2ad1.png)

そして、この四角形の画像をマスク画像と重ね、透明なところだけ NDVI の値を採用し、透明でない部分は 0 にします

```elixir
mask_mat
|> Evision.Mat.to_nx(EXLA.Backend)
|> Nx.slice_along_axis(3, 1, axis: 2)
|> Nx.tile([3])
|> Nx.select(0, Evision.Mat.to_nx(city_ndvi_bbox_map, EXLA.Backend))
|> Nx.as_type(:u8)
|> Evision.Mat.from_nx_2d()
```

![ndvi_city.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/6b81e57a-aa4f-a52d-4b57-9ecccc666871.png)

植生を藤崎町の形に切り取ることができました

## 時系列による推移

今までの一連の処理をモジュール化します

```elixir
defmodule CityNDVI do
  def get_lat_lon header_info do
    ndvi_longitudes = [
      header_info.left_top.longitude,
      header_info.right_top.longitude,
      header_info.left_bottom.longitude,
      header_info.right_bottom.longitude
    ]

    ndvi_latitudes = [
      header_info.left_top.latitude,
      header_info.right_top.latitude,
      header_info.left_bottom.latitude,
      header_info.right_bottom.latitude
    ]

    ndvi_min_longitude = Enum.min(ndvi_longitudes)
    ndvi_max_longitude = Enum.max(ndvi_longitudes)
    ndvi_min_latitude = Enum.min(ndvi_latitudes)
    ndvi_max_latitude = Enum.max(ndvi_latitudes)

    {ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude}
  end

  def wrap_rotate tensor, degree do
    {src_height, src_width} = Nx.shape(tensor)

    affine =
      {src_width / 2, src_height / 2}
      |> Evision.getRotationMatrix2D(degree, 1)
      |> Evision.Mat.to_nx(EXLA.Backend)

    cos = Nx.abs(affine[[0, 0]])
    sin = Nx.abs(affine[[0, 1]])

    dst_width =
      Nx.add(Nx.multiply(src_height, sin), Nx.multiply(src_width, cos))
      |> Nx.to_number()
      |> trunc()

    dst_height =
      Nx.add(Nx.multiply(src_height, cos), Nx.multiply(src_width, sin))
      |> Nx.to_number()
      |> trunc()

    bias =
      Nx.tensor([
        [0, 0, (dst_width - src_width) / 2],
        [0, 0, (dst_height - src_height) / 2]
      ])

    new_affine = Nx.add(affine, bias)

    Evision.warpAffine(tensor, new_affine, {dst_width, dst_height})
  end

  def correct ndvi_tensor, header_info, ndvi_lat_lon do
    {ndvi_min_longitude, ndvi_max_longitude, _, _} =
      ndvi_lat_lon

    rotated_mat = wrap_rotate(ndvi_tensor, -header_info.degree)

    diff_longitude =
      (header_info.left_top.longitude - header_info.right_top.longitude) -
      (header_info.left_bottom.longitude - header_info.right_bottom.longitude)

    {ndvi_height, ndvi_width} = rotated_mat.shape

    pix_per_longitude = ndvi_width / (ndvi_max_longitude - ndvi_min_longitude)

    diff_x = trunc(diff_longitude / 2 * pix_per_longitude)

    affine =
      [
        [1, 0, diff_x],
        [0, 1, 0]
      ]
      |> Nx.tensor()
      |> Nx.as_type(:f32)

    Evision.warpAffine(rotated_mat, affine, {ndvi_width, ndvi_height})
  end

  def get_city_bbox corrected_mat, mask_mat, city_geojson, ndvi_lat_lon do
    {ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude} =
      ndvi_lat_lon

    city_coordinates =
      city_geojson.geometries
      |> Enum.map(& &1.coordinates)
      |> List.flatten()

    city_longitudes = Enum.map(city_coordinates, &elem(&1, 0))
    city_latitudes = Enum.map(city_coordinates, &elem(&1, 1))

    city_min_longitude = Enum.min(city_longitudes)
    city_max_latitude = Enum.max(city_latitudes)

    {ndvi_height, ndvi_width} = corrected_mat.shape

    pix_per_longitude = ndvi_width / (ndvi_max_longitude - ndvi_min_longitude)
    pix_per_latitude = ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)

    {city_height, city_width, 4} = mask_mat.shape

    city_left_top_x = trunc((city_min_longitude - ndvi_min_longitude) * pix_per_longitude)
    city_left_top_y = trunc((ndvi_max_latitude - city_max_latitude) * pix_per_latitude)

    corrected_mat
    |> Evision.Mat.to_nx(EXLA.Backend)
    |> then(
      & &1[
        [
          city_left_top_y..(city_left_top_y + city_height - 1),
          city_left_top_x..(city_left_top_x + city_width - 1)
        ]
      ]
    )
  end

  def crop city_bbox, mask_mat do
    mask_mat
    |> Evision.Mat.to_nx(EXLA.Backend)
    |> Nx.slice_along_axis(3, 1, axis: 2)
    |> Nx.select(0, city_bbox)
  end

  def get_city_ndvi ndvi_tensor, header_info, mask_mat, city_geojson do
    ndvi_lat_lon = get_lat_lon(header_info)

    corrected_mat = correct(ndvi_tensor, header_info, ndvi_lat_lon)

    city_bbox = get_city_bbox(corrected_mat, mask_mat, city_geojson, ndvi_lat_lon)

    city_mean_ndvi =
      mask_mat
      |> Evision.Mat.to_nx(EXLA.Backend)
      |> Nx.slice_along_axis(3, 1, axis: 2)
      |> Nx.squeeze()
      # マスク内だけで平均する
      |> Nx.select(0, 1)
      |> then(&Nx.weighted_mean(city_bbox, &1))
      |> Nx.to_number()

    city_ndvi_map =
      city_bbox
      |> Nx.multiply(128)
      |> Nx.add(128)
      |> Nx.as_type(:u8)
      |> Evision.applyColorMap(Evision.Constant.cv_COLORMAP_WINTER())
      |> Evision.Mat.to_nx(EXLA.Backend)

    croped_city_ndvi_map =
      mask_mat
      |> Evision.Mat.to_nx(EXLA.Backend)
      |> Nx.slice_along_axis(3, 1, axis: 2)
      |> Nx.tile([3])
      |> Nx.select(0, city_ndvi_map)
      |> Nx.as_type(:u8)

    {city_mean_ndvi, croped_city_ndvi_map}
  end
end
```

各シーンについて行政区域内の NDVI を算出します

このとき、可視化のためのテンソルと、グラフ化のための平均値を併せて計算しています

```elixir
city_ndvi_list =
  ndvi_list
  |> Enum.map(fn {header_info, ndvi_tensor} ->
    # 観測日をタブにする
    date =
      header_info.datetime
      |> NaiveDateTime.to_iso8601()
      |> String.slice(0..9)

    IO.inspect(date)
    
    {city_ndvi, city_ndvi_map} =
      CityNDVI.get_city_ndvi(ndvi_tensor, header_info, mask_mat, city_geojson)

    {date, city_ndvi, city_ndvi_map}
  end)
```

日付毎の NDVI 画像をタブ表示で出力します

可視化のためのテンソルは Evision （OpenCV） の BGR 色空間になっているため、 Nx.reverse で BGR に変えてから表示しています

```elixir
city_ndvi_list
|> Enum.map(fn {date, _, city_ndvi_map} ->
  {date, city_ndvi_map |> Nx.reverse(axes: [2]) |> Kino.Image.new()}
end)
|> Kino.Layout.tabs()
```

![fujisaki_ndvi.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/2027c667-1974-3657-ca84-eee339b42ee7.gif)

冬はもちろん暗いですが、6月もまだ暗く、8月に一気に明るくなっているのが分かります

下に前回の、衛星データ全体での NDVI 推移を示します

![ndvi.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/026fdbe9-b499-0df9-96bf-b90626b6efd4.gif)

衛星データ全体で見たときは山林の範囲が広く、6月の時点で木々が青々と茂っている様子が分かります

藤崎町に限定すると田んぼの範囲が広いため、稲の育つ8月になってから NDVI が高くなっています

更にグラフで比べてみましょう

```elixir
plot_data =
  city_ndvi_list
  |> Enum.map(fn {date, city_mean_ndvi, _} ->
    %{
      date: date,
      ndvi: city_mean_ndvi
    }
  end)
```

```elixir
VegaLite.new(width: 700)
|> VegaLite.data_from_values(plot_data)
|> VegaLite.mark(:bar)
|> VegaLite.encode_field(:x, "date", type: :temporal)
|> VegaLite.encode_field(:y, "ndvi", type: :quantitative)
```

藤崎町だけに限定すると、 8 月だけ突出して NDVI が高くなっていることが分かります

![fujisaki_ndvi_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/bbdbbb8a-e949-65cc-8548-0706193e13fa.png)

下に示す衛星データ全体の NDVI では、 6 月でもかなり高くなっています

![graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/328fc767-e216-8a15-2738-23aa62317e3b.png)

## まとめ

行政区域で NDVI を切り抜くことによって、山林が多いのか、田んぼが多いのか、あるいはどちらも少ないのか、などがある程度比較できそうですね

更に田んぼ単位で切り取れば、そこが休耕地なのかどうかも分かりそうです

次回は農林水産省の筆ポリゴンデータを使って、田んぼ単位の NDVI を算出し、休耕地を探してみます

https://qiita.com/RyoWakabayashi/items/4039649daa0d22a3bda0