はじめに
経済産業省の「衛星データ利用環境整備・ソリューション開発支援事業」において、弊社は衛星データ無料利用事業者として支援を受けています
本事業では日本発の衛星データプラットフォーム「Tellus」のデータの一部を無料利用できます
（一般利用者が無料利用できるデータもあります）
「Tellus」から取得できる様々なデータを使って何をどう実装できるのか、現在調査検証中です
連載で「Tellus」の公式メディア「宙畑-sorabatake-」の記事を参考にしながら、 Elixir で衛星データを扱っています
参考にする記事
今回は植生（植物の集まり、植物がどれくらい生えているか）を見るために算出した NDVI(Normalized Difference Vegetation Index) = 正規化植生指数について、特定の市区町村内だけを切り抜きます
これによって、市区町村毎の植生の比較や、市区町村内の植生推移が可視化できます
引き続き Elixir の Livebook を使っていきます
まだ読んでいない方は環境構築編、データ取得編から読んでください
連載記事
本記事で実装したノートブックはこちら
出典
行政区域データ
「国土数値情報（行政区域データ）」（国土交通省）（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html）を加工して作成
衛星データ
提供：だいち(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 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 を取得し、画像として表示します
{header_info, ndvi_tensor} = Enum.at(ndvi_list, 0)
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})
これを地図に重ねてみます
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]
]
# 地図にプロットするために BASE64 エンコードする
get_data_url = fn mat ->
Evision.imencode(".png", mat)
|> Base.encode64()
|> then(&"data:image/png;base64,#{&1}")
end
img_url =
ndvi_map
|> Evision.resize({4960, 4960})
|> 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)
ここまでは前回のおさらいです
市区町村のポリゴンデータを用意する
国土交通省の行政区域データダウンロードページから青森県、令和4年のデータをダウンロードします
行政区域データについては以下の記事を参照してください
この記事ではダウンロードしたデータを /tmp/GML に展開したものとして進めます
gml_dir = "/tmp/GML"
geojson_file =
gml_dir
# ファイル一覧取得
|> File.ls!()
# `.geojson` で終わるもののうち先頭を取得
|> Enum.find(&String.ends_with?(&1, ".geojson"))
行政コード = 02361 (藤崎町)のデータを抽出します
行政コード一覧はコチラ
city_geojson = %Geo.GeometryCollection{
geometries: Enum.filter(geojson_data.geometries, &(&1.properties["N03_007"] == "02361"))
}
航空写真の地図に重ねてみます
city_map =
MapLibre.new(center: {140.5, 40.68}, zoom: 11, style: :terrain)
|> MapLibre.add_geo_source("city_geojson", city_geojson)
|> MapLibre.add_layer(
id: "city_geojson",
source: "city_geojson",
type: :fill,
paint: [fill_color: "#00ff00", fill_opacity: 0.5]
)
地図上の藤崎町の境界線と、緑色で表示した行政区域データが見事に重なっています
この行政区域データのポリゴン（多角形）で NDVI を切り抜き、藤崎町内の NDVI を可視化します
NDVI 画像の補正
衛星データの形状を確認してみます
{ndvi_height, ndvi_width, 3} = Evision.Mat.shape(ndvi_map)
以下のような結果が表示されます
{8000, 7244, 3}
つまり、衛星データは縦 8,000 ピクセル × 横 7,244 ピクセル × 3 色の画像データになっています
しかし、先ほど地図上に NDVI をプロットしたのを見ても分かるとおり、この画像データは地軸に対して斜めに傾いています
どれくらい傾いているかはヘッダー情報に格納されています
header_info.degree
値は以下のような度（°）で与えられます
13.849887
この分回転して地軸に対して並行にすれば、行政区域データに重ねることができそうです
四隅の角が画像内に残るように画像を回転する関数を定義します
wrap_rotate = fn mat, degree ->
# 元の大きさ
{src_height, src_width, 3} = mat.shape
# 四隅がはみ出ることを考慮しない回転のアフィン変換
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()
# (変換後サイズ - 変更前サイズ) / 2 の分だけずらして画像の中心を回転の中心に合わせる
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(mat, new_affine, {dst_width, dst_height})
end
NDVI の画像を回転させます
rotated_map = wrap_rotate.(ndvi_map, -header_info.degree)
Evision.resize(rotated_map, {640, 640})
また、画像は長方形になっていますが、緯度経度を見ると、これも少しずれていることが分かります
長方形であるならば、回転していたとしても上辺の経度の差と下辺の経度の差は同じはずです
diff_longitude =
(header_info.left_top.longitude - header_info.right_top.longitude) -
(header_info.left_bottom.longitude - header_info.right_bottom.longitude)
しかし、この上辺と下辺の経度の差を計算すると、以下のような値になります
-0.008156400000018493
この値の分だけ台形になっていそうです
正確には台形補正をすべきですが、ここでは単純に誤差の半分だけ左側に水平移動します
この辺りは正確ではないです
しかし何らかのか補正をしないとぴったり当てはまりません
衛星データに詳しい方、ご教示いただけると助かります
1ピクセルあたりの緯度経度を計算します
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}
{ndvi_height, ndvi_width, 3} = Evision.Mat.shape(rotated_map)
pix_per_longitude = ndvi_width / (ndvi_max_longitude - ndvi_min_longitude)
pix_per_latitude = ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)
{pix_per_longitude, pix_per_latitude}
計算結果は以下のようになります
{8483.693398865544, 11116.281262859115}
この値を使って、横方向への補正値を計算します
diff_x = trunc(diff_longitude / 2 * pix_per_longitude)
補正値は以下の値になります
-34
補正値を使って画像を横にずらします
affine =
[
[1, 0, diff_x],
[0, 1, 0]
]
|> Nx.tensor()
|> Nx.as_type(:f32)
corrected_map = Evision.warpAffine(rotated_map, affine, {ndvi_width, ndvi_height})
Evision.resize(corrected_map, {640, 640})
ほとんど変わっていませんが、僅かながら右上の角と右端の間に隙間ができています
その分、左下の角が画像からはみ出ました
マスク画像による切り抜き
NDVI の画像を切り抜くためのマスク画像を作成します
# 緯度経度の最大最小を求める
city_coordinates =
city_geojson.geometries
|> Enum.map(& &1.coordinates)
|> List.flatten()
city_longitudes = Enum.map(city_coordinates, &elem(&1, 0))
city_latitudes = Enum.map(city_coordinates, &elem(&1, 1))
city_min_longitude = Enum.min(city_longitudes)
city_max_longitude = Enum.max(city_longitudes)
city_min_latitude = Enum.min(city_latitudes)
city_max_latitude = Enum.max(city_latitudes)
{city_min_longitude, city_max_longitude, city_min_latitude, city_max_latitude}
# 最大最小の差からマスク画像の大きさを求める
city_width = trunc((city_max_longitude - city_min_longitude) * pix_per_longitude)
city_height = trunc((city_max_latitude - city_min_latitude) * pix_per_latitude)
{city_width, city_height}
行政区域データの各座標（緯度経度）をマスク画像内の座標（XY）に変換します
# 緯度経度の最小から最大をピクセル数の0から幅・高さにスケールする
normalized_points =
city_geojson.geometries
|> Enum.map(& &1.coordinates)
|> Enum.map(fn coordinate ->
coordinate
|> Enum.at(0)
|> Enum.map(fn {x, y} ->
[
trunc((x - city_min_longitude) * pix_per_longitude),
# 縦方向は緯度が北緯の場合逆転するため、高さから引く
trunc(city_height - (y - city_min_latitude) * pix_per_latitude)
]
end)
|> Nx.tensor(type: :s32)
end)
真っ黒で不透明な空画像を作ります
# 空画像を用意する
empty_mat =
[0, 0, 0, 255]
|> Nx.tensor(type: :u8)
|> Nx.tile([city_height, city_width, 1])
|> Evision.Mat.from_nx_2d()
Evision.resize(empty_mat, {640, 640})
空画像内に行政区域データを多角形として書き込みます
# ポリゴンを透明色で塗りつぶす
mask_mat = Evision.fillPoly(empty_mat, normalized_points, {0, 0, 0, 0})
Evision.resize(mask_mat, {640, 640})
白い部分は透明になっています
地図へのプロット
マスク画像を地図にプロットしてみます
city_center = {
(city_min_longitude + city_max_longitude) / 2,
(city_min_latitude + city_max_latitude) / 2
}
city_bbox_coordinates =
[
[city_min_longitude, city_max_latitude],
[city_max_longitude, city_max_latitude],
[city_max_longitude, city_min_latitude],
[city_min_longitude, city_min_latitude]
]
image_url =
mask_mat
|> Evision.resize({trunc(city_height / 4), trunc(city_width / 4)})
|> get_data_url.()
MapLibre.new(center: city_center, zoom: 11, style: :terrain)
|> MapLibre.add_source("city_mask", type: :image, url: image_url, coordinates: city_bbox_coordinates)
|> MapLibre.add_layer(id: "overlay", source: "city_mask", type: :raster)
藤崎町をくり抜くことができています
NDVI を行政区域で切り抜く
いよいよ NDVI を切り抜きます
まず、行政区域を含む領域を四角形で切り抜きます
# 緯度経度を元に NDVI 画像内の行政区域の位置を計算する
city_left_top_x = trunc((city_min_longitude - ndvi_min_longitude) * pix_per_longitude)
city_left_top_y = trunc((ndvi_max_latitude - city_max_latitude) * pix_per_latitude)
{city_left_top_x, city_left_top_y}
city_ndvi_bbox_map =
corrected_map
|> Evision.Mat.to_nx(EXLA.Backend)
|> then(
& &1[
[
city_left_top_y..(city_left_top_y + city_height - 1),
city_left_top_x..(city_left_top_x + city_width - 1)
]
]
)
|> Evision.Mat.from_nx_2d()
そして、この四角形の画像をマスク画像と重ね、透明なところだけ NDVI の値を採用し、透明でない部分は 0 にします
mask_mat
|> Evision.Mat.to_nx(EXLA.Backend)
|> Nx.slice_along_axis(3, 1, axis: 2)
|> Nx.tile([3])
|> Nx.select(0, Evision.Mat.to_nx(city_ndvi_bbox_map, EXLA.Backend))
|> Nx.as_type(:u8)
|> Evision.Mat.from_nx_2d()
植生を藤崎町の形に切り取ることができました
時系列による推移
今までの一連の処理をモジュール化します
defmodule CityNDVI do
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, ndvi_lat_lon do
{ndvi_min_longitude, ndvi_max_longitude, _, _} =
ndvi_lat_lon
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)
Evision.warpAffine(rotated_mat, affine, {ndvi_width, ndvi_height})
end
def get_city_bbox corrected_mat, mask_mat, city_geojson, ndvi_lat_lon do
{ndvi_min_longitude, ndvi_max_longitude, ndvi_min_latitude, ndvi_max_latitude} =
ndvi_lat_lon
city_coordinates =
city_geojson.geometries
|> Enum.map(& &1.coordinates)
|> List.flatten()
city_longitudes = Enum.map(city_coordinates, &elem(&1, 0))
city_latitudes = Enum.map(city_coordinates, &elem(&1, 1))
city_min_longitude = Enum.min(city_longitudes)
city_max_latitude = Enum.max(city_latitudes)
{ndvi_height, ndvi_width} = corrected_mat.shape
pix_per_longitude = ndvi_width / (ndvi_max_longitude - ndvi_min_longitude)
pix_per_latitude = ndvi_height / (ndvi_max_latitude - ndvi_min_latitude)
{city_height, city_width, 4} = mask_mat.shape
city_left_top_x = trunc((city_min_longitude - ndvi_min_longitude) * pix_per_longitude)
city_left_top_y = trunc((ndvi_max_latitude - city_max_latitude) * pix_per_latitude)
corrected_mat
|> Evision.Mat.to_nx(EXLA.Backend)
|> then(
& &1[
[
city_left_top_y..(city_left_top_y + city_height - 1),
city_left_top_x..(city_left_top_x + city_width - 1)
]
]
)
end
def crop city_bbox, mask_mat do
mask_mat
|> Evision.Mat.to_nx(EXLA.Backend)
|> Nx.slice_along_axis(3, 1, axis: 2)
|> Nx.select(0, city_bbox)
end
def get_city_ndvi ndvi_tensor, header_info, mask_mat, city_geojson do
ndvi_lat_lon = get_lat_lon(header_info)
corrected_mat = correct(ndvi_tensor, header_info, ndvi_lat_lon)
city_bbox = get_city_bbox(corrected_mat, mask_mat, city_geojson, ndvi_lat_lon)
city_mean_ndvi =
mask_mat
|> 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_bbox, &1))
|> Nx.to_number()
city_ndvi_map =
city_bbox
|> Nx.multiply(128)
|> Nx.add(128)
|> Nx.as_type(:u8)
|> Evision.applyColorMap(Evision.Constant.cv_COLORMAP_WINTER())
|> Evision.Mat.to_nx(EXLA.Backend)
croped_city_ndvi_map =
mask_mat
|> Evision.Mat.to_nx(EXLA.Backend)
|> Nx.slice_along_axis(3, 1, axis: 2)
|> Nx.tile([3])
|> Nx.select(0, city_ndvi_map)
|> Nx.as_type(:u8)
{city_mean_ndvi, croped_city_ndvi_map}
end
end
各シーンについて行政区域内の NDVI を算出します
このとき、可視化のためのテンソルと、グラフ化のための平均値を併せて計算しています
city_ndvi_list =
ndvi_list
|> Enum.map(fn {header_info, ndvi_tensor} ->
# 観測日をタブにする
date =
header_info.datetime
|> NaiveDateTime.to_iso8601()
|> String.slice(0..9)
IO.inspect(date)
{city_ndvi, city_ndvi_map} =
CityNDVI.get_city_ndvi(ndvi_tensor, header_info, mask_mat, city_geojson)
{date, city_ndvi, city_ndvi_map}
end)
日付毎の NDVI 画像をタブ表示で出力します
可視化のためのテンソルは Evision （OpenCV） の BGR 色空間になっているため、 Nx.reverse で BGR に変えてから表示しています
city_ndvi_list
|> Enum.map(fn {date, _, city_ndvi_map} ->
{date, city_ndvi_map |> Nx.reverse(axes: [2]) |> Kino.Image.new()}
end)
|> Kino.Layout.tabs()
冬はもちろん暗いですが、6月もまだ暗く、8月に一気に明るくなっているのが分かります
下に前回の、衛星データ全体での NDVI 推移を示します
衛星データ全体で見たときは山林の範囲が広く、6月の時点で木々が青々と茂っている様子が分かります
藤崎町に限定すると田んぼの範囲が広いため、稲の育つ8月になってから NDVI が高くなっています
更にグラフで比べてみましょう
plot_data =
city_ndvi_list
|> Enum.map(fn {date, city_mean_ndvi, _} ->
%{
date: date,
ndvi: city_mean_ndvi
}
end)
VegaLite.new(width: 700)
|> VegaLite.data_from_values(plot_data)
|> VegaLite.mark(:bar)
|> VegaLite.encode_field(:x, "date", type: :temporal)
|> VegaLite.encode_field(:y, "ndvi", type: :quantitative)
藤崎町だけに限定すると、 8 月だけ突出して NDVI が高くなっていることが分かります
下に示す衛星データ全体の NDVI では、 6 月でもかなり高くなっています
まとめ
行政区域で NDVI を切り抜くことによって、山林が多いのか、田んぼが多いのか、あるいはどちらも少ないのか、などがある程度比較できそうですね
更に田んぼ単位で切り取れば、そこが休耕地なのかどうかも分かりそうです
次回は農林水産省の筆ポリゴンデータを使って、田んぼ単位の NDVI を算出し、休耕地を探してみます