はじめに

NDVI（Normalized Difference Vegetation Index）は近赤外バンドと赤バンドの反射率から計算される指标で、植生の活性度を反映します。耕作地は生育期にNDVIが高くなり、補地や耕作放棄地は通常NDVIが低いことから、衛星画像の時系列NDVIを用いて耕作状況を分類できます。

本記事では、3時期の光学NDVIを用いて耕作地・耕作放棄地を分類するロジックと、QGISのフィールド計算機で使える計算式を紹介します。

NDVIの計算

衛星の赤バンドと近赤外バンドからNDVIを算出します。例として、sentinel2の‭2024-04-14‬のNDVIフィールドを作成する式を示します。

with_variable('red', "red_20240414" / 10000,
with_variable('nir', "nir_20240414" / 10000,
  if( ("@nir" + "@red") = 0, NULL,
      ( "@nir" - "@red" ) / ( "@nir" + "@red" ) )
))

同様に時期ごとにNDVI列を作成します。

maxNDVI・minNDVI・diffNDVIの計算

3時期のNDVIの中から最大値・最小値と、その差分（diff）を求めます。缺死を除外する安全版の式例です。列名はあなたのデータに合わせて置き換えてください。

最大値（maxNDVI）

with_variable(
  vals',
  array
  
  
  _filter(array("ndvi_20240414", "ndvi_20240802", "ndvi_20241001"), @element IS NOT NULL),
  if(array_length(@vals) = 0, NULL, maximum(@vals))
)

最小値（minNDVI）

wi
th_variable(
  'vals',
  array_filter(array("ndvi_20240414", "ndvi_20240802", "ndvi_20241001"), @element IS NOT NULL),
  if(array_length(@vals) = 0, NULL, minimum(@vals))
)

差分（diffNDVI）

if("max_ndvi" IS NULL OR "min_ndvi" IS NULL,
   NULL,
   "max_ndvi" - "min_ndvi")

ここで max_ndvi と min_ndvi は前段で計算した列名です。

分類ロジックとフィールド演算式

耕作状況の判断は次のルールとしました。

耕作放棄地：最大NDVIが 0.2 未満
耕作地：最大NDVIが 0.5 以上かつ diffNDVIが 0.25 以上
調査必要：上記以外（放棄候補や缺死を含む）

この3クラスを数値0】2に割り当てると、スタイル設定や統計処理がしやすくなります。

CASE
  WHEN "MAX" < 0.2 THEN 2        -- 耕作放棄地
  WHEN "MAX" >= 0.5 AND "DIFF" >= 0.25 THEN 0  -- 耕作地
 
  ELSE 1                         -- 調査必要
END

ここで "MAX" と "DIFF" はそれぞれ前段で計算したフィールド名です。出力フィールド型を数値型に設定しておくと分類が分かりやすくなります。

0= 耕作地、1= 調査必要、2= 耕作放棄地というラベルは必要に応じてテキストに置き換えてもかまいません。

おわりに

3時期のNDVIを使った単純な闘値分類は、光学データが収得できる地域・時期では手軽に耕作状況を把握できる手法です。一方で、雷や雷影の影響、多年生作物や作付失敗による誤判定、小区画の畑での解像度不足などに注意が必要です。

闘値は地域や作物に応じて調整する必要があります。
注：この記事はAI（ChatGPT）によって作成された作業メモです。