NDVI

livebooks/tellus/ndvi.livemd

Ryo Wakabayashi

@RyoWakabayashi

elixir-learning

More notebooks

NDVI

Mix.install(
  [
    {:nx, "~> 0.9"},
    {:evision, "~> 0.2"},
    {:exla, "~> 0.9"},
    {:kino, "~> 0.14"},
    {:kino_maplibre, "~> 0.1"},
    {:kino_vega_lite, "~> 0.1"}
  ],
  config: [nx: [default_backend: EXLA.Backend]]
)

データの指定

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"
]

衛星データを画像として表示する

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

だいち(ALOS) AVNIR-2 の仕様はこちら

https://www.eorc.jaxa.jp/ALOS/jp/alos/sensor/avnir2_j.htm

scene_id = Enum.at(scene_id_list, 0)

一つのシーンには以下のデータが含まれている

…_thumb.png: サムネイル用画像（通常のカラー画像）
…_webcog.tif: ブラウズ画像（通常のカラー画像）
HDR-…-001.txt: ヘッダーファイル
IMG-01-…-001.tif: バンド1（青色光バンド）画像
IMG-02-…-001.tif: バンド2（緑色光バンド）画像
IMG-03-…-001.tif: バンド3（赤色光バンド）画像
IMG-04-…-001.tif: バンド4（近赤外線バンド）画像

file_path_list =
  "/tmp/#{scene_id}"
  |> File.ls!()
  |> Enum.sort()
  |> Enum.map(fn filename ->
    ["/tmp", scene_id, filename]
    |> Path.join()
  end)

file_path_list
# テキストデータであるヘッダファイルは除く
|> Enum.filter(fn file_path -> Path.extname(file_path) != ".txt" end)
|> Enum.map(fn file_path ->
  file_path
  |> Evision.imread()
  # 大きすぎるのでリサイズ
  |> Evision.resize({640, 640})
end)
|> Kino.Layout.grid(columns: 2)

ヘッダーファイルの読込

ヘッダーファイルの仕様はこちら

header_path = Enum.find(file_path_list, &amp;(Path.extname(&amp;1) == ".txt"))

# ヘッダーファイルのテキストを表示してみる
File.read!(header_path)

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

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

header_info = HeaderInfo.read(header_path)

地図へのプロット

red_img_path =
  file_path_list
  |> Enum.find(fn file_path ->
    file_path
    |> Path.basename()
    # 赤色光バンドの画像
    |> String.starts_with?("IMG-03")
  end)

red_img = Evision.imread(red_img_path)

# 大きすぎるのでリサイズして表示
Evision.resize(red_img, {640, 640})

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

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

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

img_url =
  red_img
  |> Evision.resize({640, 640})
  |> 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)

カラー画像の生成

blue_img =
  file_path_list
  |> Enum.find(&amp;(&amp;1 |> Path.basename() |> String.starts_with?("IMG-01")))
  |> Evision.imread()

green_img =
  file_path_list
  |> Enum.find(&amp;(&amp;1 |> Path.basename() |> String.starts_with?("IMG-02")))
  |> Evision.imread()

[
  Evision.resize(red_img, {640, 640}),
  Evision.resize(green_img, {640, 640}),
  Evision.resize(blue_img, {640, 640})
]
|> Kino.Layout.grid(columns: 3)

bgr_img =
  [blue_img, green_img, red_img]
  |> Enum.map(fn img ->
    img
    |> Evision.Mat.to_nx(EXLA.Backend)
    # ３チャネル（同一値）になっているので、１チャネルだけ取り出す
    |> Nx.slice_along_axis(0, 1, axis: 2)
  end)
  # チャネル方向に結合
  |> Nx.concatenate(axis: 2)
  |> Evision.Mat.from_nx_2d()

Evision.resize(bgr_img, {640, 640})

img_url =
  bgr_img
  |> Evision.resize({4000, 4000})
  |> 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)

NDVIの算出

NDVI(Normalized Difference Vegetation Index) = 正規化植生指数

$ NDVI = \frac{NIR - Red}{NIR + Red} $

# 赤色光バンド
red_tensor =
  red_img
  |> Evision.Mat.to_nx(EXLA.Backend)
  # ゲインをかけてオフセットを足す
  |> Nx.multiply(header_info.red_band.gain)
  |> Nx.add(header_info.red_band.offset)
  |> Nx.slice_along_axis(0, 1, axis: 2)

# 近赤外線バンド
nir_tensor =
  file_path_list
  |> Enum.find(&amp;(&amp;1 |> Path.basename() |> String.starts_with?("IMG-04")))
  |> Evision.imread()
  |> Evision.Mat.to_nx(EXLA.Backend)
  # ゲインをかけてオフセットを足す
  |> Nx.multiply(header_info.nir_band.gain)
  |> Nx.add(header_info.nir_band.offset)
  |> Nx.slice_along_axis(0, 1, axis: 2)

ndvi_tensor =
  Nx.select(
    # 0 除算をしないため、 NIR と Red の両方が 0 でないところだけ演算する
    Nx.multiply(
      Nx.not_equal(red_tensor, 0),
      Nx.not_equal(nir_tensor, 0)
    ),
    # NDVI の演算
    Nx.divide(
      Nx.subtract(nir_tensor, red_tensor),
      Nx.add(nir_tensor, red_tensor)
    ),
    0
  )

NDVI は -1 から 1 の範囲になっているので、 0 から 255 に変換してヒートマップを表示する

ndvi_map =
  ndvi_tensor
  |> Nx.multiply(128)
  |> Nx.add(128)
  |> Nx.as_type(:u8)
  |> Evision.Mat.from_nx_2d()
  |> then(&amp;[src: &amp;1, colormap: Evision.Constant.cv_COLORMAP_WINTER()])
  |> Evision.applyColorMap()

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

img_url =
  ndvi_map
  |> Evision.resize({640, 640})
  |> 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)

時系列で確認する

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_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_list
|> Enum.map(fn {header_info, ndvi_tensor} ->
  center = {header_info.center.longitude, header_info.center.latitude}

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

  img_url =
    ndvi_tensor
    |> Nx.multiply(128)
    |> Nx.add(128)
    |> Nx.as_type(:u8)
    |> Evision.resize({640, 640})
    |> then(&amp;[src: &amp;1, colormap: Evision.Constant.cv_COLORMAP_WINTER()])
    |> Evision.applyColorMap()
    |> get_data_url.()

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

  # 観測日をタブにする
  date =
    header_info.datetime
    |> NaiveDateTime.to_iso8601()
    |> String.slice(0..9)

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

plot_data =
  ndvi_list
  |> Enum.map(fn {header_info, ndvi_tensor} ->
    avg_ndvi =
      ndvi_tensor
      |> Nx.mean()
      |> Nx.to_number()

    %{
      datetime: header_info.datetime,
      ndvi: avg_ndvi
    }
  end)

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

Back