はじめに
経済産業省の「衛星データ利用環境整備・ソリューション開発支援事業」において、弊社は衛星データ無料利用事業者として支援を受けています
本事業では日本発の衛星データプラットフォーム「Tellus」のデータの一部を無料利用できます
（一般利用者が無料利用できるデータもあります）
「Tellus」から取得できる様々なデータを使って何をどう実装できるのか、現在調査検証中です
連載で「Tellus」の公式メディア「宙畑-sorabatake-」の記事を参考にしながら、 Elixir で衛星データを扱っています
参考にする記事
今回も Elixir の Livebook を使っていくので、まだ読んでいない方は連載の環境構築編をご一読ください
連載記事
本記事で実装したノートブックはこちら
出典
行政区域データ
「国土数値情報（行政区域データ）」（国土交通省）（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v3_1.html）を加工して作成
衛星データ
提供：だいち(ALOS) AVNIR-2 データ(JAXA)
Tellus のアカウント作成
Tellus から衛星データを取得するため、まずは Tellus のアカウントを作ります
公式ガイドを参照してください
トークンの作成
Tellus の API を呼び出すためのトークンを作成します
画面右上のログインユーザー名->「開発環境」をクリックします
「トークンの発行」をクリックすると、新しいトークンが作成され、下に表示されます
これで API を呼ぶ準備ができました
衛星データの選び方（Tellus Traveler）
参考元の記事では、とりあえず AVNIR-2 のデータを使うことが前提になっていますが、そもそもどうやってデータを選べばいいのでしょうか
Tellus には Tellus Traveler という機能があり、ここでどのようなデータがあるのか探索することができます
ここでデータセットを選択したり、データの観測日や地理情報を指定することで検索が可能です
「どのデータセットか？」というのがよく分からなくても、おそらく「どこのデータか？」は決まっていると思います
例えば参考元の記事で使っている青森県藤崎町近辺を指定して検索してみましょう
すると、検索結果はこうなります
青い四角が該当する衛星データの領域なのですが、重なりすぎていて何が何だか分かりません
場所だけを指定してもデータセットを選択していないと、該当件数が多すぎて全く絞れていないのです
ではデータセットを「【Tellus公式】AVNIR-2_1B1」だけに絞ってみましょう
そうするとかなり絞り込まれますが、これでも枠が重なり過ぎていてよく分かりません
結局一件一件詳細を開かないと、使えるデータなのかどうか確かめるのは難しいです
しかもここで恣意的に選んだデータの ID をプログラムに直書きするような形になるのはあまり美しいとは言えません
というわけで、 Elixir Livebook でこの問題を解決し、明確に衛星データを選択してみます
Elixir Livebook による衛星データの選別
セットアップ
Livebook を起動し、新しいノートブック（もしくは私の実装したノートブック）を開いてください
一番先頭のセル（セットアップセル:「Notebook dependencies and setup」 と書いてある場所）に以下の内容を入力してください
キーボードの「Shift+Enter」でセルの内容を実行します
Mix.install(
[
{:nx, "~> 0.4"},
{:evision, "~> 0.1"},
{:exla, "~> 0.4"},
{:req, "~> 0.3"},
{:geo, "~> 3.4"},
{:kino, "~> 0.8"},
{:kino_maplibre, "~> 0.1.3"}
],
config: [
nx: [
default_backend: EXLA.Backend
]
]
)
ここでは以下の依存モジュールをインストールしています
- Nx: 行列演算
- Evision: 画像処理
- EXLA: 行列演算の高速化
- Req: HTTPクライアント
- Geo: 地理情報システム
- Kino: Livebook での入出力
- Kino MapLibre: 地図出力
また、 Nx による行列演算を EXLA を使って高速化する設定をしています
認証情報の設定
実行後、下のセル（黒いブロック）に以下の内容を入力し、実行します
token_input = Kino.Input.password("Token")
テキストエリアが表示されるので、 Tellus で作っておいたトークンを入力します
モジュールの作成
次のセルに以下の内容を入力します
defmodule TellusTraveler do
@base_path "https://www.tellusxdp.com/api/traveler/v1"
@data_path "#{@base_path}/datasets"
defp get_headers(token) do
%{
"Authorization" => "Bearer #{token}",
"Content-Type" => "application/json"
}
end
def get_datasets(token, is_order_required) do
url = "#{@data_path}/?is_order_required=#{is_order_required}"
headers = get_headers(token)
url
|> Req.get!(headers: headers)
|> then(& &1.body["results"])
end
def get_dataset(token, dataset_id) do
url = "#{@data_path}/#{dataset_id}/"
headers = get_headers(token)
url
|> Req.get!(headers: headers)
|> then(& &1.body)
end
def search(token, dataset_id, coordinates) do
url =
if is_list(dataset_id) do
"#{@base_path}/data-search/"
else
"#{@data_path}/#{dataset_id}/data-search/"
end
headers = get_headers(token)
request_body =
%{
intersects: %{
type: "Polygon",
coordinates: coordinates
},
query: %{},
sortby: [
%{
field: "properties.start_datetime",
direction: "asc"
}
]
}
|> Map.merge(
if is_list(dataset_id) do
%{datasets: dataset_id}
else
%{}
end
)
|> Jason.encode!()
url
|> Req.post!(body: request_body, headers: headers)
|> then(& &1.body["features"])
end
def get_data_files(token, dataset_id, data_id) do
url = "#{@data_path}/#{dataset_id}/data/#{data_id}/files/"
headers = get_headers(token)
url
|> Req.get!(headers: headers)
|> then(& &1.body["results"])
end
defp get_data_file_download_url(token, dataset_id, data_id, file_id) do
url = "#{@data_path}/#{dataset_id}/data/#{data_id}/files/#{file_id}/download-url/"
headers = get_headers(token)
url
|> Req.post!(headers: headers)
|> then(& &1.body["download_url"])
end
def download(token, dataset_id, scene_id, dist \\ "/tmp/") do
[dist, scene_id]
|> Path.join()
|> File.mkdir_p()
token
|> get_data_files(dataset_id, scene_id)
|> Enum.map(fn file ->
file_path = Path.join([dist, scene_id, file["name"]])
unless File.exists?(file_path) do
token
|> get_data_file_download_url(dataset_id, scene_id, file["id"])
|> Req.get!(output: file_path)
end
file_path
end)
end
end
Tellus Traveler の API を呼び出すためのモジュールです
API 仕様はこちら
単に API へ HTTP リクエストを投げているだけですが、ひとまず内容は気にせず実行しましょう
次に作成したモジュールを使って、データセットの一覧を取得します
お金を使いたくないので、無償利用可能なデータを取得してみます
datasets =
token_input
|> Kino.Input.read()
|> TellusTraveler.get_datasets(false)
実行すると、以下のような結果が表示されます
Tellus Traveler の API から返ってきた JSON 形式のデータです
このままでは見にくいので、必要な項目に絞ってテーブル表示します
datasets
|> Enum.map(fn dataset ->
%{
"name" => dataset["name"],
"description" => dataset["description"],
"allowed_in" => dataset["permission"]["allow_network_type"]
}
end)
|> Kino.DataTable.new(keys: ["name", "description", "allowed_in"])
34 種類、無償で使えるデータセットがありました
ちなみにこのテーブルは並び替えや検索が可能で、そのまま Excel などに貼り付けることもできます
データセットの選び方
allowed_in の項目は「どこで使えるか」を表します
- tellus: Tellus 内のみでダウンロード不可
- global: ダウンロード可能
データを加工したり分析したりしたいので、 global のデータセットを選択しましょう
また、 description の項目にデータセットの概要が書かれています
今回は植生=植物がどれくらい生えているかを知りたいので、その場合には赤色光および近赤外線のデータが必要です
これらのデータは「光学」データなので、 description に「光学」が含まれるデータセットに絞ります
opt_datasets =
datasets
|> Enum.filter(fn dataset ->
dataset["permission"]["allow_network_type"] == "global" &&
String.contains?(dataset["description"], "光学")
end)
Kino.DataTable.new(opt_datasets, keys: ["name", "description"])
これで4種類まで絞り込まれました
ちなみに、光学ではなく SAR 衛星のデータが欲しい場合は以下のようなコードになります
sar_datasets =
datasets
|> Enum.filter(fn dataset ->
dataset["permission"]["allow_network_type"] == "global" &&
String.contains?(dataset["description"], "SAR")
end)
Kino.DataTable.new(sar_datasets, keys: ["name", "description"])
シーンの選び方
シーン一覧取得
あるデータセットの、ある日時の、ある領域の衛星データを「シーン」と読んでいます
今回は青森県藤崎町の領域を含む「シーン」を選択していきます
まず、データセットは先ほど絞り込んだ無償の光学データ4種類を対象にします
dataset_id_list = Enum.map(opt_datasets, & &1["id"])
続いて、青森県藤崎町の地理情報（緯度経度で地球上のどの範囲にあるのか）を指定したいので、国土交通省の行政区域データを取得してきます
取得方法や詳細な利用方法は以下の記事を参照してください
今回は青森県の行政区域データをダウンロードし、 /tmp/GML 内に展開しておきます
データを配置したディレクトリーを指定します
gml_dir = "/tmp/GML"
データのうち、 GeoJSON ファイルを読み込みます
geojson_file =
gml_dir
# ファイル一覧取得
|> File.ls!()
# `.geojson` で終わるもののうち先頭を取得
|> Enum.find(&String.ends_with?(&1, ".geojson"))
geojson_data =
[gml_dir, geojson_file]
|> Path.join()
|> File.read!()
|> Jason.decode!()
|> Geo.JSON.decode!()
藤崎町（行政コード = 02361）だけを抜き出します
行政コードは以下のサイトを参照してください
city_geojson = %Geo.GeometryCollection{
geometries: Enum.filter(geojson_data.geometries, &(&1.properties["N03_007"] == "02361"))
}
とりあえず藤崎町を地図上にプロットしてみましょう
city_map =
MapLibre.new(center: {140.5, 40.7}, zoom: 7)
|> 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]
)
ちゃんと青森県内に藤崎町らしき形が緑色でプロットされています
Tellus Traveler API では、多角形（Polygon）で領域を指定することで、その領域を含むシーンを検索することができます
ただし、指定できる多角形の座標数は 20 までになっているため、行政区域データのようなかなり角の多い多角形は指定できません
行政区域データの緯度・経度の最大値・最小値をそれぞれ計算し、藤崎町に外接する四角形を計算しましょう
city_polygons =
city_geojson.geometries
|> Enum.at(0)
|> then(& &1.coordinates)
longitudes =
city_polygons
|> List.flatten()
|> Enum.map(&elem(&1, 0))
latitudes =
city_polygons
|> List.flatten()
|> Enum.map(&elem(&1, 1))
bbox = %{
min_longitude: Enum.min(longitudes),
max_longitude: Enum.max(longitudes),
min_latitude: Enum.min(latitudes),
max_latitude: Enum.max(latitudes),
}
city_rectangle =
[
[
[bbox.min_longitude, bbox.min_latitude],
[bbox.max_longitude, bbox.min_latitude],
[bbox.max_longitude, bbox.max_latitude],
[bbox.min_longitude, bbox.max_latitude],
[bbox.min_longitude, bbox.min_latitude],
]
]
四角形ですが、 Polygon では一周して元の位置に戻らないといけないので、最初と最後は同じ点で合計5ヶ所の緯度経度を指定しています
この四角形を使ってシーンを検索します
scenes_list =
token_input
|> Kino.Input.read()
|> TellusTraveler.search(dataset_id_list, city_rectangle)
シーン一覧が取得できました
シーン領域の可視化
シーン毎に衛星データの領域を四角形で持っています
この領域を地図上にプロットしてみましょう
scenes_list
|> Enum.map(& &1["geometry"]["coordinates"])
|> Enum.map(fn coordinates ->
coordinates
|> Enum.at(0)
|> Enum.map(& List.to_tuple(&1))
end)
|> then(fn coordinates ->
%Geo.GeometryCollection{
geometries: [%Geo.Polygon{
coordinates: coordinates
}]
}
end)
|> then(fn geojson ->
city_map
|> MapLibre.add_geo_source("area", geojson)
|> MapLibre.add_layer(
id: "area",
source: "area",
type: :line,
paint: [line_color: "#00ff00"]
)
end)
藤崎市の領域を含むシーンといっても、どうやらそれぞれ領域が異なっているようです
しかもおおよそ4種類の領域に分かれていそうです
領域によっては藤崎市にかすっているだけで、ほとんど使えなさそうなものもあるように見えます
まだまだ絞り込みが必要です
データセット毎の件数確認
データセット毎に何件取得できたか確認してみます
target_dataset_id_list =
scenes_list
# データセットIDでグルーピング
|> Enum.group_by(& &1["dataset_id"])
# それぞれの件数を取得
|> Enum.map(fn {key, value} -> {key, Enum.count(value)}end)
|> Enum.into(%{})
4種類データセットを指定していましたが、結局藤崎町の領域にデータがあるのは1種類だけでした
もしここで複数データセットがある場合は複数データセットを横断して使うのか、最もシーン数が多いデータセットを使うのかなど検討してください
データセットIDからデータセット名や概要を取得してみます
target_dataset_id_list
|> Enum.map(fn {dataset_id, count} ->
Enum.find(opt_datasets, & &1["id"] == dataset_id)
|> Map.merge(%{"count" => count})
end)
|> Kino.DataTable.new(keys: ["id", "name", "count", "description"])
というわけで選ばれたのは「【Tellus公式】AVNIR-2_1B1」でした
参考元の記事と同じものです
これ以降、このデータセットだけを対象として扱います
target_dataset_id = "ea71ef6e-9569-49fc-be16-ba98d876fb73"
データセットと雲量によるシーンの絞り込み
データセットIDと雲量を使って絞り込みます
eo:cloud_cover は 0 から 100 の範囲で雲の量を表します
この値が大きければ空が雲で覆われていて地面が観測できていない、ということになります
なので eo:cloud_cover に 25 未満を指定します
target_scenes_list =
scenes_list
|> Enum.filter(fn scene ->
scene["dataset_id"] == target_dataset_id &&
scene["properties"]["eo:cloud_cover"] < 25
end)
この時点で 31 件まで絞られます
領域によるシーンの絞り込み
続いて領域（緯度経度）で絞り込みたいと思います
同じデータセットであれば衛星で撮影できる面積は同じなので、左上座標だけを使ってグループ分けしてみましょう
Evision (OpenCV) を使って「k-平均法」によるグループ分けを行います
各シーンから左上座標を取得し、テンソルに変換します
coordinates_tensor =
target_scenes_list
|> Enum.map(fn scene ->
scene["geometry"]["coordinates"]
|> Enum.at(0)
|> Enum.at(0)
end)
|> Nx.tensor()
これを Evision.kmeans でグループ分けします
先ほど全シーンをプロットしたとき、おおよそ 4 種類の領域がありそうだったので、グループ数を 4 にしています
その他のパラメータ詳細については OpenCV の説明を参照してください
labels =
coordinates_tensor
|> Evision.kmeans(
# ４グループに分ける
4,
# 必要ないが nil 指定できないので適当なテンソルを指定する
Nx.tensor([0], type: :f32),
# 3 = TERM_CRITERIA_EPS(1) + TERM_CRITERIA_MAX_ITER(2)
{3, 10, 1.0},
# 試行回数
10,
# 中心初期化手法を使用
Evision.Constant.cv_KMEANS_PP_CENTERS
)
|> then(fn {compactness, labels, _} ->
IO.inspect("compactness: #{compactness}")
labels
|> Evision.Mat.to_nx()
|> Nx.to_flat_list()
end)
デバッグ出力している compactness が 0.01 くらいならおおよそグループ分けできています
compactness が大きよようであればグループ数を増やしてみましょう
もっと適切なグループ数を探したい場合はエルボー法を使いますが、それは後日別記事にします
何はともあれ４グループに分けることができました
出力される List の 0 1 2 3 がグループの番号になっています
このグループ番号 List を使って、シーンをグループに分けます
target_scenes_group =
[target_scenes_list, labels]
|> Enum.zip()
|> Enum.group_by(fn {_, label} -> label end, fn {scene, _} -> scene end)
グループ毎の件数を見てみます
target_scenes_group
|> Enum.map(fn {label, scenes} -> {label, Enum.count(scenes)} end)
|> Enum.into(%{})
各グループ同じくらいの件数です
これだけではピンとこないので、各グループの先頭シーンを地図上にプロットしてみます
# グループ毎に先頭シーンの領域を取得
scene_geojson_map =
target_scenes_group
|> Enum.map(fn {label, scene_list} ->
scene_list
|> Enum.at(0)
|> then(& &1["geometry"]["coordinates"])
|> Enum.at(0)
|> Enum.map(& List.to_tuple(&1))
|> then(fn coordinates ->
%Geo.GeometryCollection{
geometries: [%Geo.Polygon{
coordinates: [coordinates]
}]
}
end)
|> then(fn geojson ->
{label, geojson}
end)
end)
|> Enum.into(%{})
# 藤崎町の中心座標
center =
{
(bbox.min_longitude + bbox.max_longitude) / 2,
(bbox.min_latitude + bbox.max_latitude) / 2
}
まずシーンと藤崎町がどのように重なっているか見るため、地図上で藤崎町を黒く塗りつぶしておきます
city_map =
MapLibre.new(center: center, zoom: 7)
|> MapLibre.add_geo_source("city_geojson", city_geojson)
|> MapLibre.add_layer(
id: "city",
source: "city_geojson",
type: :fill,
paint: [fill_color: "#000000"]
)
そして、ここからがこの記事のハイライトです
各グループの領域を地図にプロットし、タブで切り替えて表示できるようにします
scene_geojson_map
|> Enum.map(fn {label, scene_geojson} ->
map =
city_map
|> MapLibre.add_geo_source("area", scene_geojson)
|> MapLibre.add_layer(
id: "area",
source: "area",
type: :fill,
paint: [fill_color: "#00ff00", fill_opacity: 0.5]
)
{Integer.to_string(label), map}
end)
|> Kino.Layout.tabs()
グループ1は藤崎市がすみっこすぎて、一部がはみ出していそうです
グループ２は藤崎市をかすっているだけで使えそうにありません
グループ0が一番良さそうです
こんな可視化ができるのが Livebook のすごいところですね
Kino.Layout についてはこちら
では、グループ内の各シーンが本当に藤崎町を含んでいるのかも確認してみましょう
target_scenes_group[0]
|> Enum.map(fn scene ->
coordinates =
scene["geometry"]["coordinates"]
|> Enum.at(0)
|> Enum.map(& List.to_tuple(&1))
scene_geojson =
%Geo.GeometryCollection{
geometries: [%Geo.Polygon{
coordinates: [coordinates]
}]
}
map =
city_map
|> MapLibre.add_geo_source("area", scene_geojson)
|> MapLibre.add_layer(
id: "area",
source: "area",
type: :fill,
paint: [fill_color: "#00ff00", fill_opacity: 0.5]
)
# 観測日をタブにする
date = String.slice(scene["properties"]["start_datetime"], 0..9)
{date, map}
end)
|> Kino.Layout.tabs()
どのシーンもほぼ同じ領域を示していますね
このグループのデータを使えば、同じ領域の推移を観測できそうです
シーンIDの一覧を見てみましょう
scene_id_list =
target_scenes_group[0]
|> Enum.map(& &1["id"])
以下の結果が表示されます
["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"]
このシーンIDは参考元記事で使ったノートブックで指定しているものとほぼ同じ（一つだけこちらが多い）です
したがって、 Elixir Livebook を使うことで、 Tellus Traveler の画面検索と同等（もしくはより根拠を示せる形で）衛星データを選択できる、ということになります
衛星データのダウンロード
選択したシーンIDを使ってデータをダウンロードしましょう
データが大きいので、相当の時間が掛かるのと、ストレージが満杯にならないように気をつけてください
scene_id_list
|> Enum.map(fn scene_id ->
token_input
|> Kino.Input.read()
|> TellusTraveler.download(target_dataset_id, scene_id)
end)
ダウンロードした画像の表示
scene_id = Enum.at(scene_id_list, 2)
"/tmp/#{scene_id}"
|> File.ls!()
|> Enum.filter(fn filename -> Path.extname(filename) != ".txt" end)
|> Enum.sort()
|> Enum.map(fn filename ->
["/tmp", scene_id, filename]
|> Path.join()
|> Evision.imread()
# 大きすぎるのでリサイズ
|> Evision.resize({640, 640})
end)
|> Kino.Layout.grid(columns: 2)
十和田湖や陸奥湾が見えますね
下4枚の画像は各波長帯の光学データで、左下が赤色光、右下が赤外線のデータです
次回以降、これらの衛星データを使って植生を計算し、地図へのプロットや季節による推移のグラフ化を実装していきます
まとめ
Elixir Livebook を使うことで、かなり視覚的、効果的に衛星データの選別を行うことができました
特にタブ表示はかなり強力です
今後、 Livebook がより広く使われることを期待します
次回は NDVI という値を計算することによって、植物の分布（植生）を可視化します