衛星データプラットフォーム 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/d65dd3a508284e67c8cd)
- 休耕地探索編 <- ここ

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

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

## 出典

### 筆ポリゴンデータ

「筆ポリゴンデータ（2022年度公開）」（農林水産省）（<https://www.maff.go.jp/j/tokei/porigon/>）を加工して作成

### 衛星データ

提供：だいち(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 算出用のモジュールを用意します

[NDVI算出編](https://qiita.com/RyoWakabayashi/items/b4fa2f7efa2a5d249a44)で使ったモジュールに補正処理やヒートマップ化まで加えています

```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 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) do
    {ndvi_min_longitude, ndvi_max_longitude, _, _} = get_lat_lon(header_info)

    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)

    rotated_mat
    |> Evision.warpAffine(affine, {ndvi_width, ndvi_height})
    |> 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
      )
      |> correct(header_info)

    {header_info, ndvi_tensor}
  end

  def to_heatmap(ndvi_tensor) do
    ndvi_tensor
    |> Nx.multiply(128)
    |> Nx.add(128)
    |> Nx.as_type(:u8)
    |> Evision.applyColorMap(Evision.Constant.cv_COLORMAP_WINTER())
  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_tensor
|> NDVI.to_heatmap()
|> Evision.resize({640, 640})
```

![corrected_ndvi.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/fdb771d2-75aa-4d54-5cb1-b23ad3f556bb.png)

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

## 筆ポリゴンデータを用意する

[農林水産省の筆ポリゴン公開サイト](https://open.fude.maff.go.jp/)から青森県藤崎町の2022年度のデータをダウンロードします

筆ポリゴンデータについては以下の記事を参照してください

https://qiita.com/RyoWakabayashi/items/4d5ff1a2ae8644d3a7b2

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

```elixir
json_file = "/tmp/2022_023612.json"
```

```elixir
geojson_data =
  json_file
  |> File.read!()
  |> Jason.decode!()
  |> Geo.JSON.decode!()
```

航空写真の地図に重ねてみます

```elixir
# 地図に表示する中心座標を求める
longitudes =
  geojson_data.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()
  |> Enum.map(&elem(&1, 0))

latitudes =
  geojson_data.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()
  |> Enum.map(&elem(&1, 1))

center = {
  (Enum.min(longitudes) + Enum.max(longitudes)) / 2,
  (Enum.min(latitudes) + Enum.max(latitudes)) / 2
}
```

まず全ての農地を表示してみます

```elixir
city_map =
  MapLibre.new(center: center, zoom: 11, style: :terrain)
  |> MapLibre.add_geo_source("geojson", geojson_data)
  |> MapLibre.add_layer(id: "line", source: "geojson", type: :line)
```

![fujisaki_all_fields_map.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/ec5a4b2b-cf99-4c23-2213-53c0e9dbb3f4.gif)

藤崎町内の田んぼが綺麗に線で区切られています

今回は田んぼを対象にするので、絞り込みます

```elixir
fields_geojson_data =
  geojson_data.geometries
  |> Enum.filter(&(&1.properties["land_type"] == 100)) # 田んぼを指定
  |> then(fn geometries ->
    %Geo.GeometryCollection{geometries: geometries}
  end)

fields_map =
  MapLibre.new(center: center, zoom: 11, style: :terrain)
  |> MapLibre.add_geo_source("geojson", fields_geojson_data)
  |> MapLibre.add_layer(id: "line", source: "geojson", type: :line)
```

![fujisaki_rice_map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/50a67072-a7aa-1dd8-8c5c-f13594976259.png)

## NDVI を田んぼの形に切り抜く

まず、全田んぼの形で NDVI を切り抜いてみます

```elixir
# 全田んぼを含む領域の緯度経度の最大最小
fields_coordinates =
  fields_geojson_data.geometries
  |> Enum.map(& &1.coordinates)
  |> List.flatten()

fields_longitudes = Enum.map(fields_coordinates, &elem(&1, 0))
fields_latitudes = Enum.map(fields_coordinates, &elem(&1, 1))

fields_min_longitude = Enum.min(fields_longitudes)
fields_max_longitude = Enum.max(fields_longitudes)
fields_min_latitude = Enum.min(fields_latitudes)
fields_max_latitude = Enum.max(fields_latitudes)

{fields_min_longitude, fields_max_longitude, fields_min_latitude, fields_max_latitude}
```

```elixir
# NDVI 画像のサイズ
{ndvi_height, ndvi_width} = Nx.shape(ndvi_tensor)
```

```elixir
# NDVI 画像の緯度経度の最大最小
{ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude} =
  NDVI.get_lat_lon(header_info)
```

```elixir
# 緯度経度あたりのピクセル数
{pix_per_longitude, pix_per_latitude} =
  {
    ndvi_width / (ndvi_max_longitude - ndvi_min_longitude),
    ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)
  }
```

```elixir
# 最大最小の差から全田んぼを含む領域のサイズを求める
fields_width = trunc((fields_max_longitude - fields_min_longitude) * pix_per_longitude)
fields_height = trunc((fields_max_latitude - fields_min_latitude) * pix_per_latitude)

{fields_width, fields_height}
```

```elixir
# NDVI 画像内での相対的な位置を求める
fields_lft_top_x = trunc((fields_min_longitude - ndvi_min_longitude) * pix_per_longitude)
fields_lft_top_y = trunc((ndvi_max_latitude - fields_max_latitude) * pix_per_latitude)

{fields_lft_top_x, fields_lft_top_y}
```

```elixir
# NDVI 画像から全田んぼを含む領域を切り抜く
city_tensor =
  ndvi_tensor[
    [
      fields_lft_top_y..(fields_lft_top_y + fields_height - 1),
      fields_lft_top_x..(fields_lft_top_x + fields_width - 1)
    ]
  ]

city_tensor
|> NDVI.to_heatmap()
```

![fields_area_ndvi.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/5f9589ab-b696-1c45-ebec-a555541581ac.png)

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

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

```elilir
# NDVI 画像にポリゴンを重ねる
empty_mat
|> Evision.fillPoly(normalized_points, {1})
|> Evision.Mat.to_nx(EXLA.Backend)
|> Nx.slice_along_axis(3, 1, axis: 2)
|> Nx.tile([3])
|> Nx.select(0, city_tensor |> NDVI.to_heatmap() |> Evision.Mat.to_nx(EXLA.Backend))
|> Nx.as_type(:u8)
|> Evision.Mat.from_nx_2d()
```

![fields_masked_ndvi.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/16770b71-8f25-f692-48aa-97c3a857ad6d.png)

NDVI 画像から田んぼの領域だけを切り抜きました

## 田んぼ毎の平均NDVI算出

田んぼ毎に平均 NDVI を算出したいので、田んぼ毎のポリゴン内部の NDVI だけを対象に平均を計算します

```elixir
fileds_mean_ndvi_tensor =
  normalized_points
  |> Enum.map(fn field_points ->
    empty_mat
    |> Evision.fillPoly([field_points], {1})
    |> 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_tensor, &1))
    |> Nx.new_axis(0)
  end)
  |> Nx.concatenate()
```

上記のコードの実行結果は以下のようになります

```elixir
#Nx.Tensor<
  f32[6463]
  EXLA.Backend<host:0, 0.1947203998.3918659595.62433>
  [-0.030122317373752594, -5.798943457193673e-4, 0.11528561264276505, -0.03583264350891113, -0.035946596413850784, -0.05206907168030739, 0.0724615529179573, 0.0020273893605917692, 3.7007121136412024e-4, 0.2709633708000183, 0.19496579468250275, 0.13714201748371124, -0.028841190040111542, 0.036340270191431046, -0.02724451571702957, 0.04559248313307762, -0.029057519510388374, -0.03282160684466362, 0.0018153662094846368, -0.03287540376186371, 0.032594528049230576, 0.08571335673332214, -0.044919345527887344, -0.029049117118120193, -0.044206585735082626, -0.038011737167835236, 0.18794618546962738, -0.04767170548439026, -0.017481816932559013, 0.03601334989070892, 0.1939351111650467, -0.0317654088139534, -0.013227726332843304, -0.029473701491951942, 0.009326216764748096, 0.07990346848964691, 0.1634574830532074, -0.010438849218189716, -0.036008432507514954, -0.04567558318376541, -0.0075868344865739346, -0.05032256618142128, -0.05560816079378128, -0.00643538311123848, 0.07092132419347763, -0.005250267684459686, -0.04875883460044861, 0.2746807038784027, 0.009388083592057228, 0.10280899703502655, ...]
>
```

田んぼの土地が全部で 6463 筆あり、それぞれの領域内の平均 NDVI が算出できました

NDVI は -1 から 1 の間になりますが、ぱっと見小さい値ばかりのようです

これは先頭シーンが6月で、まだお米が育っていないからだと思われます

数字ではいまいち分かりづらいので、田んぼ毎にNDVI画像の当該領域を表示してみましょう

NDVI 値の高い田んぼと低い田んぼをそれぞれ表示します

```elixir
city_map =
  city_tensor
  |> NDVI.to_heatmap()
  |> Evision.Mat.to_nx(EXLA.Backend)

# NDVI の高い順にインデックスを取得
ordered = Nx.argsort(fileds_mean_ndvi_tensor, direction: :desc)

# 上位５つと下位5つを取得
top_and_worst =
  [ordered[0..4], ordered[-5..-1]]
  |> Nx.concatenate()
  |> Nx.to_flat_list()

top_and_worst
|> Enum.map(fn index ->
  field_points = Enum.at(normalized_points, index)

  [min_x, min_y] = field_points |> Nx.reduce_min(axes: [0]) |> Nx.to_flat_list()
  [max_x, max_y] = field_points |> Nx.reduce_max(axes: [0]) |> Nx.to_flat_list()

  empty_mat
  |> Evision.fillPoly([field_points], {1})
  |> Evision.Mat.to_nx(EXLA.Backend)
  |> Nx.slice_along_axis(3, 1, axis: 2)
  |> Nx.tile([3])
  |> Nx.select(0, city_map)
  |> then(& &1[[min_y..max_y, min_x..max_x]])
  |> Nx.as_type(:u8)
  |> Evision.Mat.from_nx_2d()
  |> Evision.resize({(max_x - min_x + 1) * 100, (max_y - min_y + 1) * 100})
  |> Evision.Mat.to_nx(EXLA.Backend)
  |> Nx.reverse(axes: [2])
  |> Kino.Image.new()
end)
|> Kino.Layout.grid(columns: 5)
```

![ndvi_top_worst.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/8c736a29-7c5d-ccd9-46f9-b29daf5e6961.png)

上段が NDVI 上位５つの田んぼで、下段が NDVI 下位5つの田んぼです

今回使っているヒートマップでは値が高いほど緑、低いほど青にしているため、確かに高いところと低いところが表示できています

## 時系列による推移

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

```elixir
defmodule FieldNDVI do
  def get_fields_lat_lon(fields_geojson_data) do
    # 緯度経度の最大最小を求める
    fields_coordinates =
      fields_geojson_data.geometries
      |> Enum.map(& &1.coordinates)
      |> List.flatten()

    fields_longitudes = Enum.map(fields_coordinates, &elem(&1, 0))
    fields_latitudes = Enum.map(fields_coordinates, &elem(&1, 1))

    {
      Enum.min(fields_longitudes),
      Enum.max(fields_longitudes),
      Enum.min(fields_latitudes),
      Enum.max(fields_latitudes)
    }
  end

  def get_pix_per_lat_lon(ndvi_tensor, header_info) do
    {ndvi_height, ndvi_width} = Nx.shape(ndvi_tensor)

    {ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude} =
      NDVI.get_lat_lon(header_info)

    {
      ndvi_width / (ndvi_max_longitude - ndvi_min_longitude),
      ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)
    }
  end

  def get_city_tensor(ndvi_tensor, header_info, fields_geojson_data) do
    {
      fields_min_longitude, fields_max_longitude,
      fields_min_latitude, fields_max_latitude
    } =
      get_fields_lat_lon(fields_geojson_data)

    {ndvi_min_longitude, _, _, ndvi_max_latitude} =
      NDVI.get_lat_lon(header_info)

    {pix_per_longitude, pix_per_latitude} =
      get_pix_per_lat_lon(ndvi_tensor, header_info)

    {fields_width, fields_height} =
      {
        trunc((fields_max_longitude - fields_min_longitude) * pix_per_longitude),
        trunc((fields_max_latitude - fields_min_latitude) * pix_per_latitude)
      }

    {fields_lft_top_x, fields_lft_top_y} = 
      {
        trunc((fields_min_longitude - ndvi_min_longitude) * pix_per_longitude),
        trunc((ndvi_max_latitude - fields_max_latitude) * pix_per_latitude)
      }

    ndvi_tensor[
      [
        fields_lft_top_y..(fields_lft_top_y + fields_height - 1),
        fields_lft_top_x..(fields_lft_top_x + fields_width - 1)
      ]
    ]
  end

  def get_normalized_points(ndvi_tensor, header_info, fields_geojson_data) do
    {
      fields_min_longitude, fields_max_longitude,
      fields_min_latitude, _
    } =
      get_fields_lat_lon(fields_geojson_data)

    {pix_per_longitude, pix_per_latitude} =
      get_pix_per_lat_lon(ndvi_tensor, header_info)

    fields_height =
      trunc((fields_max_longitude - fields_min_longitude) * pix_per_longitude)

    fields_geojson_data.geometries
    |> Enum.map(& &1.coordinates)
    |> Enum.map(fn coordinate ->
      coordinate
      |> Enum.at(0)
      |> Enum.map(fn {x, y} ->
        [
          trunc((x - fields_min_longitude) * pix_per_longitude),
          # 縦方向は緯度が北緯の場合逆転するため、高さから引く
          trunc(fields_height - (y - fields_min_latitude) * pix_per_latitude)
        ]
      end)
      |> Nx.tensor(type: :s32)
    end)
  end

  def get_field_mean_ndvi(field_points, city_tensor, empty_mat) do
    empty_mat
    |> Evision.fillPoly([field_points], {1})
    |> 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_tensor, &1))
    |> Nx.new_axis(0)
  end

  def get_fields_mean_ndvi(ndvi_tensor, header_info, fields_geojson_data) do
    normalized_points =
      get_normalized_points(ndvi_tensor, header_info, fields_geojson_data)
    
    city_tensor =
      get_city_tensor(ndvi_tensor, header_info, fields_geojson_data)

    {fields_height, fields_width} = Nx.shape(city_tensor)

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

    normalized_points
    |> Enum.map(&get_field_mean_ndvi(&1, city_tensor, empty_mat))
    |> Nx.concatenate()
  end
end
```

全シーンの全田んぼについて平均NDVIを算出します

```elixir
fields_mean_ndvi_tensor =
  ndvi_list
  |> Enum.map(fn {header_info, ndvi_tensor} ->
    IO.inspect(header_info.datetime)
    FieldNDVI.get_fields_mean_ndvi(ndvi_tensor, header_info, fields_geojson_data)
  end)
  |> Nx.stack()
```

実行結果は以下のようなテンソルになります

```elixir
#Nx.Tensor<
  f32[8][6463]
  EXLA.Backend<host:0, 0.1947203998.3918921740.173733>
  [
    [0.01952052302658558, -0.0012143654748797417, 0.10302069783210754, -0.017521314322948456, -0.0379776731133461, -0.058998268097639084, 0.07677389681339264, -1.6859121387824416e-4, 0.013595733791589737, 0.2902262210845947, 0.25921204686164856, 0.19830869138240814, 0.007779088336974382, 4.1169748874381185e-4, -0.046043623238801956, 0.020258240401744843, -0.03633059561252594, -0.05142595246434212, 5.380002039601095e-5, -0.019140347838401794, 0.02145955339074135, 0.03939266875386238, -0.03585371747612953, -0.014490923844277859, -0.0888376235961914, -0.010234634391963482, -0.017845263704657555, -0.04945629462599754, -0.02470400370657444, 2.236171712866053e-4, 0.12306370586156845, -0.005346057470887899, -0.03845560923218727, -0.03238437697291374, 0.039897769689559937, 0.05320005118846893, 0.04643165320158005, -0.003063499927520752, -0.02353859320282936, -0.03602850064635277, 0.001951472251676023, -0.02061539888381958, -0.04269783943891525, 0.009131945669651031, 0.04980583116412163, 0.03606180474162102, -0.05278889834880829, 0.1688486486673355, 0.0673358216881752, -0.004329879768192768, ...],
    ...
  ]
>
```

8 がシーン数、 6463 が田んぼの数です

このテンソルを先頭の次元で区切れば、各田んぼの時系列変化を見ることができます

今回の目的は休耕地を探すことです

耕地であれば冬にNDVIが低く、夏にNDVIが高くなるはずですが、休耕地の場合は植物がないため、一貫してNDVIが低くなります

つまり、NDVIの時間による変化の大きさを見れば休耕地を探索できそうです（最大値を見る、という手もありますが）

どれくらい値にばらつきがあるか、というのは標準偏差で表すことができます

$ \bar{x} $ を $ x_1 $ から $ x_n $ の平均としたとき、標準偏差は以下のような式で表します

$ \sigma = \sqrt{\frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})^2} $

つまり、標準偏差とは平均と各値との差を積み重ねたものです

平均から離れた値が多ければ標準偏差が大きくなり、どの値も平均に近ければ標準偏差は小さくなります

Elixir では `Nx.standard_deviation` に `axes: [0]` を指定することで、先頭の次元方向（時系列）に対して標準偏差を求めることができます

```elixir
fields_std_ndvi_tensor = Nx.standard_deviation(fields_mean_ndvi_tensor, axes: [0])
```

実行結果は以下のような、田んぼの数の1次元テンソルになります

```elixir
#Nx.Tensor<
  f32[6463]
  EXLA.Backend<host:0, 0.1947203998.3918921740.173747>
  [0.12534421682357788, 0.2096170037984848, 0.1888715773820877, 0.20814335346221924, 0.24141794443130493, 0.19366270303726196, 0.18087920546531677, 0.26057735085487366, 0.16534438729286194, 0.08822018653154373, 0.14118514955043793, 0.1646028608083725, 0.14609003067016602, 0.1719924509525299, 0.26297447085380554, 0.057613153010606766, 0.2437814474105835, 0.2538285553455353, 0.23053482174873352, 0.22712987661361694, 0.23607677221298218, 0.22217749059200287, 0.23310844600200653, 0.1816251575946808, 0.17127683758735657, 0.24108855426311493, 0.23602735996246338, 0.11905814707279205, 0.24868209660053253, 0.22345130145549774, 0.13466420769691467, 0.20992381870746613, 0.23128917813301086, 0.2358301430940628, 0.09428323805332184, 0.14384020864963531, 0.24990765750408173, 0.2318825125694275, 0.23812319338321686, 0.23407939076423645, 0.21120834350585938, 0.21793213486671448, 0.23806947469711304, 0.2251603901386261, 0.21261337399482727, 0.20978417992591858, 0.25278839468955994, 0.17720484733581543, 0.22474639117717743, 0.2297637164592743, ...]
>
```

標準偏差の大きい順に並べた上位５つ、下位５つを見比べてみましょう

```elixir
ordered = Nx.argsort(fields_std_ndvi_tensor, direction: :desc)
```

各田んぼのNDVIの推移をグラフ化してみます

```elixir
plot_ndvi = fn mean_ndvi_tensor, title ->
  plot_data =
    mean_ndvi_tensor
    |> Nx.to_flat_list()
    |> Enum.map(fn index ->
      ndvi_list
      |> Enum.with_index()
      |> Enum.map(fn {{header_info, _}, scene_index} ->
        %{
          datetime: header_info.datetime,
          index: index,
          ndvi: Nx.to_number(fields_mean_ndvi_tensor[[scene_index, index]])
        }
      end)
    end)
    |> List.flatten()

  VegaLite.new(width: 700, title: title)
  |> VegaLite.data_from_values(plot_data)
  |> VegaLite.mark(:line)
  |> VegaLite.encode_field(:x, "datetime", type: :temporal)
  |> VegaLite.encode_field(:y, "ndvi", type: :quantitative, scale: %{"domain" => [-1.0, 1.0]})
  |> VegaLite.encode_field(:color, "index", type: :nominal)
end
```

```Elixir
[
  {ordered[[0..4]], "top"},
  {ordered[[-5..-1]], "worst"}
]
|> Enum.map(fn {tensor, title} -> plot_ndvi.(tensor, title) end)
|> Kino.Layout.grid()
```

上位５つの田んぼはどれも8月に NDVI が高くなり、それ以外は低くなっています

![top_ndvi_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/20f9af85-78f8-bff3-542a-85c89b22d79e.png)

下位５つの田んぼは常に NDVI が 0　付近です

![worst_ndvi_graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/9ff93c07-ff90-c486-69ae-88471bebdf90.png)

今度は田んぼ毎の NDVI 推移を画像で確認してみましょう

```elixir
display_field_ndvi = fn index_list, ndvi_tensor, header_info, fields_geojson_data ->
  normalized_points =
    FieldNDVI.get_normalized_points(ndvi_tensor, header_info, fields_geojson_data)
  
  city_tensor =
    FieldNDVI.get_city_tensor(ndvi_tensor, header_info, fields_geojson_data)

  city_map =
    city_tensor
    |> NDVI.to_heatmap()
    |> Evision.Mat.to_nx(EXLA.Backend)

  {fields_height, fields_width} = Nx.shape(city_tensor)

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

  index_list
  |> Enum.map(fn index ->
    field_points = Enum.at(normalized_points, index)

    [min_x, min_y] = field_points |> Nx.reduce_min(axes: [0]) |> Nx.to_flat_list()
    [max_x, max_y] = field_points |> Nx.reduce_max(axes: [0]) |> Nx.to_flat_list()

    empty_mat
    |> Evision.fillPoly([field_points], {1})
    |> Evision.Mat.to_nx(EXLA.Backend)
    |> Nx.slice_along_axis(3, 1, axis: 2)
    |> Nx.tile([3])
    |> Nx.select(0, city_map)
    |> then(& &1[[min_y..max_y, min_x..max_x]])
    |> Nx.as_type(:u8)
    |> Evision.Mat.from_nx_2d()
    |> Evision.resize({(max_x - min_x + 1) * 100, (max_y - min_y + 1) * 100})
    |> Evision.Mat.to_nx(EXLA.Backend)
    |> Nx.reverse(axes: [2])
    |> Kino.Image.new()
  end)
  |> Kino.Layout.grid(columns: 5)
end
```

```elixir
top_and_worst =
  [ordered[0..4], ordered[-5..-1]]
  |> Nx.concatenate()
  |> Nx.to_flat_list()

ndvi_list
|> Enum.map(fn {header_info, ndvi_tensor} ->  
  date =
    header_info.datetime
    |> NaiveDateTime.to_iso8601()
    |> String.slice(0..9)

  display = display_field_ndvi.(top_and_worst, ndvi_tensor, header_info, fields_geojson_data)

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

![ndvi_fields_map.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1485835/3d178288-6905-96d1-263d-f0f6dfb72daf.gif)

上段が上位５つ、下段が下位５つです

上段は8月で緑色になるのに対して、下段は常に青いままです

## まとめ

田んぼ毎に NDVI の時間推移が明確に違っていました

NDVI の標準偏差を使うことで休耕地を探すことができそうです

今回で一旦連載は終了ですが、今後も宙畑の記事を ELixir 化していきます