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

はじめに

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

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

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

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

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

参考にする記事

今回は植生（植物の集まり、植物がどれくらい生えているか）を見るため、 NDVI(Normalized Difference Vegetation Index) = 正規化植生指数を算出します

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

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

連載記事

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

出典

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

前提条件

データ取得編で対象としたシーンのデータを Tellus からダウンロードし、 /tmp/<シーンID>/ 配下に保存しているものとします

セットアップ

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

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

Mix.install(
  [
    {:nx, "~> 0.4"},
    {:evision, "~> 0.1"},
    {:exla, "~> 0.4"},
    {:kino, "~> 0.8"},
    {:kino_maplibre, "~> 0.1.3"},
    {:kino_vega_lite, "~> 0.1.7"}
  ],
  config: [
    nx: [
      default_backend: EXLA.Backend
    ]
  ]
)

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

Nx: 行列演算
Evision: 画像処理
EXLA: 行列演算の高速化
Kino: Livebook での入出力
Kino MapLibre: 地図出力
Kino VegaLite: グラフ出力

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

データの指定

データ取得編で選択したシーンIDリストを指定しておきます

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

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

前回の最後にやりましたが、改めて衛星データを画像として表示してみます

まず先頭のシーンを表示対象にします

scene_id = Enum.at(scene_id_list, 0)

対象データの保存先ディレクトリーからファイル一覧を取得します

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

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

["/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/ALAV2A040842780-OORIRFU-D069P3-20061031-001_thumb.png",
 "/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/ALAV2A040842780_webcog.tif",
 "/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/HDR-ALAV2A040842780-OORIRFU-D069P3-20061031-001.txt",
 "/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/IMG-01-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif",
 "/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/IMG-02-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif",
 "/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/IMG-03-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif",
 "/tmp/202ce08d-ba4b-4ffe-8165-109fd3a8b917/IMG-04-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif"]

ファイル名を見ただけでは何のデータか分からないので、説明書を確認してみましょう

ALOS AVNIR-2 オルソ補正画像プロダクトプロダクトフォーマット説明書 Version 2.0

これを見ると、格納しているファイルは以下のものだと書かれています

ヘッダファイル
オルソデータファイル（バンド1）
オルソデータファイル（バンド2）
オルソデータファイル（バンド3）
オルソデータファイル（バンド4）
ブラウズ画像

ヘッダファイルはテキストで衛星データの座標情報などを持っています

ブラウズ画像は Tellus Traveler 上で表示する際のカラー画像のようです

ブラウズ画像は TIFF と PNG の2種類が入っています

ではオルソデータファイル（バンド1〜4）はどういう意味でしょうか

ALOS　AVNIR-2　の製品ページを見てみましょう

すると、観測波長帯の項目に以下の記述があります

Band1：0.42 ～ 0.50µm
Band2：0.52 ～ 0.60µm
Band3：0.61 ～ 0.69µm
Band4：0.76 ～ 0.89µm

どうやらバンド1〜4はそれぞれの波長帯で観測したデータのようです

波長帯とは、、、？という方は以下の宙畑の記事を参照してください

要するに、以下のように対応しています

バンド1: 青色光
バンド2: 緑色光
バンド3: 赤色光
バンド4: 近赤外線

このうち、 NDVI の算出には赤色光と近赤外線のデータを使用します

では、これらのデータを画像としてみてみましょう

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)

各画像は以下のものを示します

左上: ブラウズ画像(PNG)
右上: ブラウズ画像(TIFF)
左中: 青色光データ
右中: 緑色光データ
左下: 赤色光データ
右下: 近赤外線データ

ヘッダーファイルの読込

この画像だけだと、どこのいつのデータなのか分からないため（実はTIFF画像に情報が埋め込まれていますが）、ヘッダーファイルの情報を読み込みます

header_path = Enum.find(file_path_list, & Path.extname(&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)

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

%HeaderInfo{
  blue_band: %BandInfo{gain: 0.588, offset: 0.0},
  green_band: %BandInfo{gain: 0.573, offset: 0.0},
  red_band: %BandInfo{gain: 0.502, offset: 0.0},
  nir_band: %BandInfo{gain: 0.557, offset: 0.0},
  center: %Coordinate{latitude: 40.5403552, longitude: 140.7569056},
  left_top: %Coordinate{latitude: 40.9672806, longitude: 140.4512007},
  right_top: %Coordinate{latitude: 40.81202, longitude: 141.2864133},
  left_bottom: %Coordinate{latitude: 40.2662829, longitude: 140.231684},
  right_bottom: %Coordinate{latitude: 40.1125883, longitude: 141.0587402},
  affine: %Affine{a: -23.9378926, b: 97.0926223, c: 64502.5142872, d: 451205.2954717},
  conversion: %Conversion{x: 10.0, y: 10.0},
  degree: 13.849887,
  map_degree: 0.0027576,
  datetime: ~N[2006-10-31 01:34:11],
  product_id: "OORIRFU"
}

データの座標や観測日時が分かりますね

地図へのプロット

では取得したヘッダー情報を元に衛星データを地図上にプロットしてみましょう

まず赤色光データを画像として読み込みます

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})

赤色光のデータですが、 BGR (Evision = OpenCV では RGB ではなく BGR です) の全てに同じ値が入っているため、灰色になって表示されます

これを地図に表示するために BASE64 エンコードします

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

MapLibre を使って、地図上に画像を重ねて表示します

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)

表示された地図を見ると、灰色の画像（衛星データ）が斜めになって青森県の航空写真に重なっています

地図をダブルクリックしてズームしていくと、衛星データと航空写真の海岸線がピッタリ重なっています

ヘッダー情報の座標を指定したことで、ピッタリ位置も縮尺も合いました

カラー画像の生成

NDVI とは関係ありませんが、せっかくなので青と緑も合わせてカラー画像を生成してみましょう

まず赤、緑、青のデータを並べてみます

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

green_img =
  file_path_list
  |> Enum.find(& &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 の順に重ねればカラー画像になります

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)

ズームしていくと、航空写真では緑色になっている田んぼが、衛星データでは茶色になっていることが分かります

そう、ヘッダー情報にある通り、このデータは 2006-10-31 に観測したものです

10月末では既にお米は収穫済で、田んぼは土が露出しているのです

NDVIの算出

いよいよ NDVI を算出します

NDVI の計算式は以下のとおりです

$ NDVI = \frac{NIR - Red}{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(& &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 を算出します

この計算は行列演算なので、各画素（ピクセル）に対して実行していますが、衛星データの左右両端には真っ黒な部分（値が 0 になっている部分）があります

値が 0 のピクセルについてこの計算をすると分母が 0 になってしまうため、その場合は計算結果を 0 とするようにします

Nx.select は第1引数の条件式に従って、真の場合は第2引数、偽の場合は第3引数の値を各要素に返してくれます

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 の値は NIR と Red の値によって以下のように変化します

NIR=255(最大),Red=0(最小) の場合: 1(最大)
NIR=0(最小),Red=255(最大) の場合: -1(最小)

従って、 NDVI の各画素は -1 〜 1 の範囲に分布しています

これを画像として表示する場合、 -1 〜 1 を 0 〜 255 に置き換えないといけません

また、灰色だと視覚的に分かりづらいため、カラーマップ（画素値と色の対応表）を 0 のとき青、 255 のとき緑となる COLORMAP_WINTER に指定します

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})

厳密にいうと　NDVI = 1 のときに値が 255 を超えてしまうのですが、ここでは無視します

十和田湖と陸奥湾の水面は真っ青になっています

また、山は緑で町は青いです

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

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)

青森県の11月末の植生の様子が可視化できました

時系列で確認する

11月末のデータを見ただけでは、本当にこれで植生が見えているのか判断できません

植物の多さを表すのであれば、春夏に高く（緑色）、秋冬に低く（青色）になるはずです

というわけで、他のシーンについても NDVI を計算し、時系列で確認してみましょう

今までやってきた NDVI の計算をモジュール化します

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 の計算を実行します

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})
    |> Evision.applyColorMap(Evision.Constant.cv_COLORMAP_WINTER())
    |> 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)

こちらでも夏が高くて冬が低いことが分かりますね

まとめ

衛星データから植生を可視化することができました

時系列の変化も非常に分かりやすいですね

次回は NDVI を町の形に切り取り、町内の植生を見てみます