はじめに
経済産業省の「衛星データ利用環境整備・ソリューション開発支援事業」において、弊社は衛星データ無料利用事業者として支援を受けています
本事業では日本発の衛星データプラットフォーム「Tellus」のデータの一部を無料利用できます
（一般利用者が無料利用できるデータもあります）
「Tellus」から取得できる様々なデータを使って何をどう実装できるのか、現在調査検証中です
連載で「Tellus」の公式メディア「宙畑-sorabatake-」の記事を参考にしながら、 Elixir で衛星データを扱っています
参考にする記事
今回は植生（植物の集まり、植物がどれくらい生えているか）を見るために算出した NDVI(Normalized Difference Vegetation Index) = 正規化植生指数について、田んぼ毎に時系列で値を算出します
これによって、田んぼでお米が育てられているのかどうか、つまり休耕地なのかどうかが判定できます
引き続き Elixir の Livebook を使っていきます
まだ読んでいない方は環境構築編、データ取得編から読んでください
連載記事
本記事で実装したノートブックはこちら
出典
筆ポリゴンデータ
「筆ポリゴンデータ（2022年度公開）」（農林水産省）（https://www.maff.go.jp/j/tokei/porigon/）を加工して作成
衛星データ
提供：だいち(ALOS) AVNIR-2 データ(JAXA)
前提条件
データ取得編 で対象としたシーンのデータを Tellus からダウンロードし、 /tmp/<シーンID>/ 配下に保存しているものとします
セットアップ
Livebook を起動し、新しいノートブック（もしくは私の実装したノートブック）を開いてください
セットアップセルに以下の内容を入力して実行してください
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: 行列演算
- Evision: 画像処理
- EXLA: 行列演算の高速化
- Geo: 地理情報システム
- Jason: JSONパーサー
- 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"
]
ヘッダーファイルの読込モジュール
NDVI算出編で使ったヘッダーファイルの読込モジュールを用意します
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
NDVIの算出モジュール
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 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 を算出します
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 を取得し、画像として表示します
{header_info, ndvi_tensor} = Enum.at(ndvi_list, 0)
ndvi_tensor
|> NDVI.to_heatmap()
|> Evision.resize({640, 640})
ここまでは前回のおさらいです
筆ポリゴンデータを用意する
農林水産省の筆ポリゴン公開サイトから青森県藤崎町の2022年度のデータをダウンロードします
筆ポリゴンデータについては以下の記事を参照してください
この記事ではダウンロードしたデータを /tmp/2022_023612.json に展開したものとして進めます
json_file = "/tmp/2022_023612.json"
geojson_data =
json_file
|> File.read!()
|> Jason.decode!()
|> Geo.JSON.decode!()
航空写真の地図に重ねてみます
# 地図に表示する中心座標を求める
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
}
まず全ての農地を表示してみます
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)
藤崎町内の田んぼが綺麗に線で区切られています
今回は田んぼを対象にするので、絞り込みます
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)
NDVI を田んぼの形に切り抜く
まず、全田んぼの形で NDVI を切り抜いてみます
# 全田んぼを含む領域の緯度経度の最大最小
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}
# NDVI 画像のサイズ
{ndvi_height, ndvi_width} = Nx.shape(ndvi_tensor)
# NDVI 画像の緯度経度の最大最小
{ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude} =
NDVI.get_lat_lon(header_info)
# 緯度経度あたりのピクセル数
{pix_per_longitude, pix_per_latitude} =
{
ndvi_width / (ndvi_max_longitude - ndvi_min_longitude),
ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)
}
# 最大最小の差から全田んぼを含む領域のサイズを求める
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}
# 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}
# 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()
# 田んぼの各ポリゴンについて、緯度経度の最小から最大をピクセル数の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)
# 空画像を用意する
empty_mat =
[0, 0, 0, 255]
|> Nx.tensor(type: :u8)
|> Nx.tile([fields_height, fields_width, 1])
|> Evision.Mat.from_nx_2d()
# 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()
NDVI 画像から田んぼの領域だけを切り抜きました
田んぼ毎の平均NDVI算出
田んぼ毎に平均 NDVI を算出したいので、田んぼ毎のポリゴン内部の NDVI だけを対象に平均を計算します
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()
上記のコードの実行結果は以下のようになります
#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 値の高い田んぼと低い田んぼをそれぞれ表示します
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 上位５つの田んぼで、下段が NDVI 下位5つの田んぼです
今回使っているヒートマップでは値が高いほど緑、低いほど青にしているため、確かに高いところと低いところが表示できています
時系列による推移
今までの一連の処理をモジュール化します
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を算出します
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()
実行結果は以下のようなテンソルになります
#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] を指定することで、先頭の次元方向（時系列）に対して標準偏差を求めることができます
fields_std_ndvi_tensor = Nx.standard_deviation(fields_mean_ndvi_tensor, axes: [0])
実行結果は以下のような、田んぼの数の1次元テンソルになります
#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, ...]
>
標準偏差の大きい順に並べた上位５つ、下位５つを見比べてみましょう
ordered = Nx.argsort(fields_std_ndvi_tensor, direction: :desc)
各田んぼのNDVIの推移をグラフ化してみます
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
[
{ordered[[0..4]], "top"},
{ordered[[-5..-1]], "worst"}
]
|> Enum.map(fn {tensor, title} -> plot_ndvi.(tensor, title) end)
|> Kino.Layout.grid()
上位５つの田んぼはどれも8月に NDVI が高くなり、それ以外は低くなっています
下位５つの田んぼは常に NDVI が 0 付近です
今度は田んぼ毎の NDVI 推移を画像で確認してみましょう
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
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()
上段が上位５つ、下段が下位５つです
上段は8月で緑色になるのに対して、下段は常に青いままです
まとめ
田んぼ毎に NDVI の時間推移が明確に違っていました
NDVI の標準偏差を使うことで休耕地を探すことができそうです
今回で一旦連載は終了ですが、今後も宙畑の記事を ELixir 化していきます