Jupyter NotebookでTellusを使ってみた〜雪質解析(1)解析準備編～
Tellusに搭載されているAPIを引っ張ることができるようになったので、雪質の解析をしながらJupyter Notebookの機能をご紹介していきたいと思います！
記事作成時から、Tellusからデータを検索・取得するAPIが変更になっております。該当箇所のコードについては、以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/access/traveler_api_20220310_
firstpart.html
2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html
本記事はTellus API公開前にリリースしたため、現在の公開APIに対応しておりません。
APIの取得方法など最新情報はTellusの開発環境でAPI引っ張ってみたをご覧ください。
寒いですね。毎年この時期になると、一旦PCを閉じてスノボに出かけたくなります。毎年スノボに行っていると、年によって雪がもはや降っていない・カチカチ・大雪など、色々な雪質のスキー場に出会います。
できる限り雪質の良いところでスノボがしたいので、衛星データで雪質を推定できるようになったら、予定もうまく立てられそう。
というわけで、雪質を推定する未来のために！Tellusを利用し、衛星データから雪を判別していきたいと思います。
衛星データプラットフォームであるTellusは、今回の最新リリースで、Pythonを使ったデータ分析で便利なJupyter Notebookに対応しました。
Tellusに搭載されているAPIを引っ張ることができるようになったので、雪質の解析をしながらJupyter Notebookの機能をご紹介していきたいと思います！
※本記事の内容を実践するためにはTellusの利用登録が必要になります。利用登録はこちらから
※本記事に掲載中のAPI呼び出しコードは過去のものになります。最新の呼び出し方法についてはTellus公式のAPIリファレンスをご参照ください。
https://www.tellusxdp.com/ja/dev/api
Tellusには、大きく2種類のアクセス方法があります。
1つはブラウザ上に表示した地図上に様々なデータを重ね、操作できるTellus OS、もう1つはプログラミングでデータを解析する開発環境(以下、Jupyter Notebook)です。
Tellusの開発環境では、各種データとJupyter NotebookがAPIで連携しています。そのため、Jupyter Notebook上からプログラミングによりAPIで引っ張ることで、データを引っ張り、解析することができます。
(1) Jupyter Notebookを起動してみる
Jupyter NotebookをTellusで使用できるようにするための手続きはここで確認できます。
https://www.tellusxdp.com/ja/tutorial
まず、Tellusで会員登録します。
「お問い合わせ導線から問い合わせせよ」とのことなので、お問い合わせ画面でカテゴリ「開発環境についての問い合わせ」を選択し、お問い合わせ内容欄に「Jupyter Notebook環境を使用したいです。」と記載し送信します。
少し待つと、窓口の方からメールが届きます。Tellusのアカウント情報を求められたので、先ほど登録したアカウント名を送信します。
その後、Jupyter Notebookのログイン情報が届きますので、初回パスワードを設定して開発環境にログインします。メールに、その後のログインの仕方が記載されているので、指示通りに操作し、ログインします。
Jupyter Notebookにログインした状態のままTellus OSにログインし、左上部のターミナルアイコンをクリックすると「Jupyter Notebookを開く」というウィンドウが現れるのでクリックします。
そうすると、Jupyter Notebookトップに遷移します。
(2) Jupyter Notebookのデフォルトディレクトリを見てみる
Jupyter Notebook環境にログインすると、いくつかのディレクトリがあります。それぞれのディレクトリは下記のような構成になっています。
●/examples：編集不可のサンプルコードが入っているようです。今回の雪質解析は、ここで試すことができます！
●/tellus_data：Tellus OSからJupyter Notebookに送信したデータが置かれるディレクトリです。Tellus OSでの送信方法はこちらの記事で読むことができます。tiff形式で送信されるため、Jupyter Notebookのプレビューでは今の所表示できないようです。
●/work ：好きなnotebookを作ることができるディレクトリです。自分の好きなnotebookを持ってくることもできます。
work以外のディレクトリは読み取り専用モードになっているので、何か新しく作りたかったり、保存したかったりするものは全てwork配下に作成することになります。
また、基本的なライブラリはプリセットでインストールされていますが、右のNewからTerminalを開けば、pip、condaコマンドでライブラリの追加インストールができます。
※但し、work、tellus_data配下以外は、アップデートなどでリセットされる可能性があるので、一旦手元に入れたものをメモしておいたほうが良さそうです。
(3) どうやって解析するか？
つぎに、Tellusに搭載される衛星データのうち、使えそうなデータを探します。
データカタログ｜Tellus
様々なスキー場を将来的に見たいので、光学画像は全国のマップがあるAVNIR-2(ALOS)、Landsatあたりが良さそうです。光学画像だけだと、雪が降っている＝雲がかかっている時の雪は見られないこともあるので、マイクロ波で雲を通過して画像を取得できる合成開口レーダ（SAR）画像もあると良さそうです。SAR画像としては、PALSAR-1(ALOS)やPALSAR-2(ALOS-2)が全国を撮像しているので、使いたいところです。
現時点のTellusでAPIとして公開されているのは、AVNIR-2だけのようなので、まずはAVNIR-2を使って解析をしてみたいと思います。※TellusのJupyter Notebook環境では、サンプルデータとして一部のPALSAR-2の衛星画像も引っ張ることができるようになっています。
(4) APIで画像取得
それでは、早速これを使ってAPIで画像取得をしてみます。APIは下記がこれから更新されていくようです。
APIリファレンス｜Tellus
というわけで、富士山周辺の光学画像をTellusから取得したいと思います。
AVNIR-2のデータは、陸域観測衛星だいち（ALOS）に搭載されていた光学センサです。青・緑・赤の光の波長帯に含まれる画像を合成することで、人の目で見るのと同じような可視光画像を作ることができます。
波長についてはこちらの記事で詳しく解説しています。
光の波長って何？ なぜ人工衛星は人間の目に見えないものが見えるのか
さて、今回用意した雪質解析のサンプルを開くと、API取得から解析までの流れを実行しながら辿ることができます。ですので、一緒に開きながら読んでいただけると理解が深まるかもしれません。
examplesディレクトリを開くと下記のような構成になっているので、まずはREADME.mdを必ず読みましょう。notebooksディレクトリには、光学、SARそれぞれの雪質解析のNotebookを記載しています。dataには、今回の解析に使用しているデータを格納しています。
README.md：今回のサンプルの全体構成などを説明しています
/notebooks：光学、SARそれぞれの雪質解析のPythonコードを、Jupyter Notebook形式で格納しています
/data：今回の解析に使用しているデータを格納しています
下記のように位置を指定して引っ張ることができます。（※サンプルコードでは、内部APIになっています。）
# APIのドメイン
URL_DOMAIN = "https://gisapi.opf-dev.jp/"
# 富士山の位置
Z = 13
X = 7252
Y = 3234
def get_data(img_type, domain=URL_DOMAIN, z=Z, x=X, y=Y):
return io.imread("{}/{}/{}/{}/{}.png".format(domain, img_type, z, x, y))
# ベースマップとなるOpen Street Map
img_osm = get_data("osm")
# 衛星データ
img_band1= get_data("av2ori/band1")
img_band2 = get_data("av2ori/band2")
img_band3 = get_data("av2ori/band3")
img_band4 = get_data("av2ori/band4")
(5) いざ解析！
以下のサンプルコードでは、画像を引っ張って、様々な基礎分析をしています。
ここでは、光学画像の基礎的な分析として、雪のセグメンテーションを取り扱います。
注意：
このnotebookで利用するTellus APIには一部サンプル目的のものが含まれます。これは予告なく変更される場合があります。
本記事での解析手順は、学術的知見に基づいたわけではなく、画像としてどのような変化が見えそうかを確認した結果になりますので正確性は保証しかねますのでご了承ください。
import numpy as np
from skimage import io
from skimage import color
from skimage import img_as_ubyte
from skimage import filters
import matplotlib.pyplot as plt
%matplotlib inlineTellusから光学画像の取得と前処理
富士山周辺の光学画像をTellusから取得します。今回はタイル化されたPNG形式のデータを扱います。位置の指定はXYZ方式を採用しています。詳細は こちら（外部サイト） をご覧ください。
以下ではALOSという衛星のAVNIR-2と呼ばれるセンサから取得されたデータを扱っています。AVNIR-2については こちら（外部サイト） をご覧ください。
注意:
ここで利用しているAPIはサンプル目的のものです。Tellus Jupyter環境のみから利用できます。正式なAPIは今後追加される予定です。
# APIのドメイン
URL_DOMAIN = "https://gisapi.opf-dev.jp/"
# 富士山の位置
Z = 13
X = 7252
Y = 3234
def get_data(img_type, domain=URL_DOMAIN, z=Z, x=X, y=Y):
return io.imread("{}/{}/{}/{}/{}.png".format(domain, img_type, z, x, y))
img_osm = get_data("osm")
img_band1= get_data("av2ori/band1")
img_band2 = get_data("av2ori/band2")
img_band3 = get_data("av2ori/band3")
img_band4 = get_data("av2ori/band4")
OpenStreetMapの確認
衛星画像だけでは、どの場所を参照しているのかわかりにくいため、参照用としてOpenStreetMapの情報も引き出せるようにします。
print(img_osm.shape)
io.imshow(img_osm)
(256, 256, 3)
上記のコードを実行すると、下記の結果が返ってきます。
AVNIR-2の確認
AVNIR-2では4つの異なる波長のデータを利用できます。大まかに、Band1は青の波長、Band2は緑の波長、Band3は赤の波長、Band4は近赤外の波長に対応します。それぞれ単独で見るためのコードは以下のようになります。
io.imshow(np.hstack((img_band1,img_band2,img_band3,img_band4)))
実行すると以下のようになります。
それぞれ可視光の青・緑・赤・（近赤外）に対応しているため、RGB（赤緑青）として1枚の画像に合成すると人の目で見た色に近いものになります。このようにして合成されたものtrue color画像と呼びます。
true colorのRGB合成
●R: Band3（赤）
●G: Band2（緑）
●B: Band1（青）
true colorの画像を出力するためのコードは以下になります。
img_true = np.c_[img_band3[:,:,0:1], img_band2[:,:,1:2], img_band1[:,:,2:3]]
io.imshow(img_true)
実行すると下記のようになります。
true color以外にも合成方法はあり、例えば植生域だけを際立たせたい時は、植物の分布域が緑で表現されるnatural colorと呼ばれる合成が用いられます。これは、近赤外線は植物の反射率が高いことを利用し、RGBの、赤色として赤の波長であるBand3、緑色として近赤外の波長であるBand4、青色として緑の波長であるBand2を割り当てたものです。
natural colorのRGB合成
●R: Band3（赤）
●G: Band4（近赤外）
●B: Band2（緑）
natural colorを合成するためのコードは以下のようになります。
img_natural = np.c_[img_band3[:,:,0:1], img_band4[:,:,0:1], img_band2[:,:,1:2]]
io.imshow(img_natural)
こちらを実行すると、以下のようになります。
直感的に白っぽいところを雪として判別してみる
では、雪を検出してみましょう。
ここでは、山の白っぽいところは雪なのでは？という直感に従って雪の抽出を行ってみます。true color画像を2つの方法でグレースケール化、更に閾値を決めて雪が抽出できるか検証してみましょう。
グレースケール化による雪判定
グレースケール化①
skimage.color.rgb2gray を使って直接グレースケール化を試みます。コードは以下のようになります。
# カラー画像からGrayscale画像への変換
img_gray_01 = color.rgb2gray(img_true)
# 値のレンジを変更 ([0, 1] -> [0, 255])
img_gray = img_as_ubyte(img_gray_01)
print("変換前: [0, 1]")
print(img_gray_01)
print("変換後: [0, 255]")
print(img_gray)
io.imshow(img_gray)
変換前: [0, 1]
[[ 0.10348667 0.10629216 0.09759294 ..., 0.13229098 0.4027102
0.28838863]
[ 0.10993098 0.10432 0.10712549 ..., 0.12278118 0.33378078
0.18417294]
[ 0.09675961 0.10320392 0.10684275 ..., 0.11662745 0.33016471
0.21639373]
...,
[ 0.10853922 0.10825647 0.1090898 ..., 1. 1. 1. ]
[ 0.10714039 0.10825647 0.1090898 ..., 1. 1. 1. ]
[ 0.11189529 0.11441804 0.11189529 ..., 1. 1. 1. ]]
変換後: [0, 255]
[[ 26 27 25 ..., 34 103 74]
[ 28 27 27 ..., 31 85 47]
[ 25 26 27 ..., 30 84 55]
...,
[ 28 28 28 ..., 255 255 255]
[ 27 28 28 ..., 255 255 255]
[ 29 29 29 ..., 255 255 255]]上記を実行すると、以下の画像を得られます。確かにグレースケール化されています。
グレースケール化②
別のグレースケール化方法も試してみましょう。一度RGB色空間からYIQ色空間へ変換します。YIQ形式は、グレースケール情報がカラーデータから分離しているため、同じ信号をカラーと白黒の両方で使用可能です。
（※ グレースケール化のアルゴリズムによっては img_gray_01 と img_yiq[:, :, 0] が等しくなりますが、skimageでは異なります）
コードにすると以下のようになります。
img_yiq = color.rgb2yiq(img_true)
img_conb = np.concatenate(
(img_yiq[:, :, 0], img_yiq[:, :, 1], img_yiq[:, :, 2]), axis=1)
io.imshow(img_conb)
実行結果は以下になります。
2つの手法で出力したグレースケール画像を比較してみましょう。比較するコードは以下になります。
# skimage.color.rgb2gray と比較
img_conb2 = np.concatenate((img_yiq[:, :, 0], img_gray_01), axis=1)
io.imshow(img_conb2)
比較した結果は以下のようになります。
雪質を目視で確認しやすくするために、白黒を反転した画像を表示します。
# 反転画像も確認
img_nega = 255 - img_gray
io.imshow(img_nega)
グレースケール化されたデータの統計情報の確認
雪のあるところを算出するために、ヒストグラムで分布を確認してみます。
print('pixel sum', np.sum(img_gray[:, :]))
print('pixel mean', np.mean(img_gray[:,:]))
print('pixel variance', np.var(img_gray[:,:]))
print('pixel stddev', np.std(img_gray[:,:]))
pixel sum 5346313
pixel mean 81.5782623291
pixel variance 5034.10730885
pixel stddev 70.9514433176ヒストグラムの確認
hists, bins = np.histogram(img_gray, 255, [0, 255])
plt.plot(hists)
出力した結果は以下のようになります。
明らかなピークは3箇所で 25 と 55 と 250 付近と分かります。 25〜55 は暗めなので森や岩、250~255 はほぼ白なので雪と推定することができます。なお、今回は目視で雲がないことを確認できているため、雪と推定することが可能ですが、雲がある場合には気を付ける必要があります。
閾値を決めて雪と思われる部分を抽出
谷の一番深い位置と思われる 200 を閾値として2値化を行います。
その処理は以下のようになります。
# 閾値処理
height, width= img_gray.shape
img_binary_gs200 = np.zeros((height, width)) # 0で初期化
img_binary_gs200 = img_gray > 200 # 閾値にもとづいて2値化
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_gs200*255), img_true), axis=1))
算出結果は以下のようになります。
一番明るいところは抽出できているものの、日陰の雪と思われる部分は全く抽出できていないようです。閾値を 80 に変更して再度2値化を行います。
img_binary_gs80 = np.zeros((height, width)) # 0で初期化
img_binary_gs80 = img_gray > 80 # 閾値にもとづいて2値化
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_gs80*255), img_true), axis=1))
雪と思われる部分をそれなりに抽出できたように見えますが、まだ右上の部分が抽出できていないため、再度閾値を 40 に下げて2値化を行ってみます。
img_binary_gs40 = np.zeros((height, width)) # 0で初期化
img_binary_gs40 = img_gray > 40 # 閾値にもとづいて2値化
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_gs40*255), img_true), axis=1))
Out[15]:
今度は範囲を過剰に取ってしまったようです。
なお、Filterプラグインを用いることで、自動で閾値を設定して2値化することもできるので、このように閾値の調整が必要な場合には非常に有用です。
th = filters.threshold_otsu(img_gray) # 閾値を自動で判別
img_binary_gsauto = np.zeros((height, width))
img_binary_gsauto = img_gray > th
print("閾値：", th)
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_gsauto*255), img_true), axis=1))
閾値： 138
Out[16]:
HSV画像による雪質判定
グレースケール画像から、雪がありそうな箇所をなんとなくは抽出できましたが、その閾値設定は非常に難しいことがわかります。
そこで、雪の検出を行う上で、グレースケール以外の手法も確認してみましょう次はHSV空間（色相・彩度・明度）に変換した画像を利用して雪の範囲を推定してみます。
# RGB空間からHSV空間へ変換
img_hsv = color.rgb2hsv(img_true)
io.imshow(img_hsv)
HSV（色相・彩度・明度）をそれぞれ確認
雪があることが、どの因子に寄与するか確認するために、色相・彩度・明度をそれぞれ確認してみましょう。
im_conb_hsv = np.concatenate(( img_hsv[:, :, 0], img_hsv[:, :, 1], img_hsv[:,:,2]), axis=1)
io.imshow(im_conb_hsv)
明度を掘り下げてみます。ヒストグラムで明度の分布を確認します。
hists, bins=np.histogram(img_hsv[:,:,2]*255, 255, [0, 255])
plt.plot(hists)
color.rgb2hsv の利用に起因すると思われる周期が見られます。ピークを分けるよう閾値を 150 に設定し、2値画像を見てみます。
img_v = img_hsv[:,:,2]
# 閾値処理
height, width= img_v.shape
img_binary_v150 = np.zeros((height, width)) # 0で初期化
img_binary_v150 = img_v > 150 / 255 # 閾値にもとづいて2値化
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_v150*255), img_true), axis=1))
Out[20]:
100 のピークの左側に閾値を設定すべく、閾値を 90 にして再度2値化を行ってみます。
img_v = img_hsv[:,:,2]
# 閾値処理
height, width= img_v.shape
img_binary_v90 = np.zeros((height, width)) # 0で初期化
img_binary_v90 = img_v > 90 / 255 # 閾値にもとづいて2値化
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_v90*255), img_true), axis=1))
RGBグレースケール化に基づく手法と明度に基づく手法の比較
さて、大きく2つの方法で雪があるかどうかの算出を試みてみましたが、それぞれの手法で算出した結果を比較してみましょう。
io.imshow(np.concatenate(
(color.gray2rgb(img_binary_gs80*255),
color.gray2rgb(img_binary_v90*255),
img_true), axis=1))
Out[22]:
明度に基づく手法はグレースケール化に基づく手法と比較して、日陰になった部分も含めて雪判定ができたようです。凹凸のある山間部の状況を確認するうえで、日陰になる箇所を考慮することは大切ですね。
更に異なる手法も考えられますが、雪判定はここまでとします。
(6) まとめ
全3回を予定しているTellus利用解析連載の雪質解析編。
第一回はJupyter Notebookにログインして、雪部分の抽出を行いました。
今回は科学的知見や理論式に基づいた解析したわけではありませんのでデータの正確性は保証できかねますが、今回行った実験的な画像解析からであっても、おおよその雪質を判定できることが分かったかと思います。
宙畑では、このほかにもTellusを用いて様々な実験や検証を引き続き行っていきたいと思います！
面白いアイディアは積極的に挑戦していこうと思いますので、分析してほしいことがありましたら、お気軽にお問い合わせいただければと思います。
皆様もこのTellusを利用して、様々な情報調べてみてはいかがでしょうか！
今回は可視光だけでしたが、次回はSARなども含めてやっていきたいと思います！
次回は積雪比率の推定に用いられることがある、NDSI(積雪比率の推定)という積雪指標を用いて雪質を分類してみたいと思います。
果たして今回の指標よりも雪質が分かりやすくなるのか？お楽しみに。
Tellusの利用登録はこちらから。
