ウェーブレット多重解像度解析(MRA)をSAR画像に適応してみた(python使用)
前回までのあらすじ
前回は、mueller行列に平均化処理をかけてノイズを除去することによって、zoomの背景画面にも使えるような合成開口レーダ(以下、SAR)画像の可視化に成功しました。
きちんと前処理したSAR画像はとても綺麗ですが、いくつか弱点があります。例えば、以下のようなことが考えられてしまいます。
・比較的平らな土地ってどこに分布してるの？(あるいは湖とか)
・写真が諏訪湖周辺なので森林が多いこともあるけど、民家がある市街地とかってどこらへんにあるのかいまいち分からない(緑一色だから森に見えちゃう)
今回は上記の問題を解決するために、SAR画像にウェーブレット多重解像度解析(MRA)という処理をかけて低周波成分と高周波成分を取り出して、平らな土地や湖、および住宅地の分布を可視化したいと思います。
提案手法としては、以下の論文を参考にして実装しました。ミュラー平均化行列を施して抽出したSAR画像は前回の記事で生成しているので、今回はその時のSAR画像(諏訪湖周辺)をそのまま使用します(実は上下反転のままなんですが笑)。今回は諏訪湖の写真も表示します。
今までのあらすじ
ちなみに、今まで取り組んできたSAR画像の生データの取り出し方や散乱特性の可視化方法などについては、以下の過去の記事を参考にしてください。
そもそも多重解像度解析とは？
多重解像度解析については真面目に解説するとそれだけでブログが一記事書けてしまうので、詳しい説明は以下を参考にしていただければと思います。超簡単に言ってしまうと、異なる解像度(周波数帯域)に分解して、SAR画像の場合は散乱特性の性質を同時に可視化できるということです。 主な使用用途としては、エッジ強調やノイズ除去に利用されることが多いです。
実際のソースコード
SAR画像のフルポラリメトリ解析に基づいて、全ての偏波成分(HH偏波、HV偏波、VV偏波)に関して多重解像度解析を行い、それぞれの結果をRGB形式で可視化することにしました。
解析に利用したソースコードは以下の通りです。
from PIL import Image import matplotlib.pyplot as plt import numpy as np import pandas as pd from PIL import ImageOps # ウェーブレット変換を行うためにPyWaveletモジュールを使いました import pywt def value(data): data = 10*np.log10(data.real**2 + data.imag**2) return data def wave_plot(HH_data,HV_data,VV_data): data = np.stack([HH_data,HV_data,VV_data],axis=2) amax = data.max() amin = data.min() #scale = 255.0/(amax-amin) #img = data*scale # ウェーブレット変換の時はスケーリング処理がいらなかったです img = data img = img.astype(np.uint8) img = Image.fromarray(img).convert('RGB') print(np.max(img),np.min(img)) plt.figure(figsize=(10,10)) img.save('wavelet.jpg', quality=95) plt.imshow(img) imgHH = value(Shh) imgHV = value(Shv) imgVV = value(Svv) wp = pywt.WaveletPacket2D(data=imgHH, wavelet='db1', mode='symmetric') hh = wp['aaa'].data # MRAを3回施した時のLL(低周波成分)を取得しています # MRAを1回施した時のHH(高周波成分)をみたい場合はhh = wp['d'].dataで取得できます wp = pywt.WaveletPacket2D(data=imgHV, wavelet='db1', mode='symmetric') hv = wp['aaa'].data wp = pywt.WaveletPacket2D(data=imgVV, wavelet='db1', mode='symmetric') vv = wp['aaa'].data wave_plot(hh,hv,vv) print(np.max(hh),np.min(hh))
結果
元画像はこちらになります(再喝)
では、実際に多重解像度解析を行なった結果と比較してみます。
MRAを3回施した時のLL(低周波成分)
三段階の多重解像度解析を施した後の低周波成分を取り出した結果がこちらになります。赤色の部分がHH偏波に対して多重解像度解析を施した低周波成分の結果になっていて、主に平坦な土地や湖などが抽出できていることが確認できました。
元々、HH偏波のデータは比較的に地表被覆物による受信電力の差が明確です。この受信電力のみを用いて画像分類を用いてもそれなりの精度が確認できるほどなので、今回は低周波成分を取り出すことで、はっきりとした地表被覆の認識が可能になりました。
MRAを1回施した時のHH(高周波成分)
こちらはかなり見辛いのですが、よく見ると画像内の白いもやもやの部分が見えるかと思います。これが住宅地のような建物が確認できる地域を示していて、諏訪湖の周辺などを白いもやが覆っていることからそれが確認できるかと思います。少し見辛い結果になってしまいましたが、多重解像度解析によって建物のある地域の特定に成功しました。
まとめ
今回はSAR画像に多重解像度解析を施すことによって、平坦な土地や建物が確認できる地域の可視化に成功しました。SAR画像はいろんな信号処理や画像処理の手法を適用できるので、いろんな解析を試してみたい方にはおすすめなので、興味がある方はぜひTellusにアカウントを登録して、専用の開発環境で衛星画像の解析をしていただけたらと思います！
実は先月から本当に色々なことがありまして、収入の事情で学生起業して個別事業主になったり、(全く別の分野で)一般科学誌向けにジャーナルを書いていたり、新規の研究のリサーチ(これもブログで触れてきた分野と全く異なる畑)もしていて、それに拍車をかけて請け負ってる仕事もあるので全然書く暇がないのですが、なんとか時間を作って解析をして行きたいと思いますのでよろしくお願いします(起業などに関してのエピソードは、また別記事でお話できたらと考えております)
次は差分干渉をするか、本格的に機械学習を試すか、またはSAR画像以外のデータを試して統合してみるか...汗 まだ未定ですが、何かしら面白い解析ができたらと思います。
PALSAR2を使った多偏波SARの解析手法(ミュラー行列+マルチルック処理)
はじめに
お久しぶりです、新星 竹です。久しぶりにSAR解析をしました。 取り組んでいる研究の解析やらジャーナルの執筆やらで佳境に入っていて、尚且つ新しい委託のお仕事を探していたり決まった後にフルコミットしたりしていた関係で、SAR画像解析をする時間が中々作れませんでした汗
今までは、TellusのGPUコンピューティング環境のセットアップからSAR画像のバイナリデータを解析してRGBで散乱特性の可視化などを行なってきました。今回は、SAR画像が解析に使えるように前処理工程を実装します。せっかくフルポラリメトリ解析ができる環境がTellusにはあるので、今回は電波が入射した表面の偏波散乱特性を表現するミュラー行列をSAR画像のデータから算出して、ミュラー行列の特徴としてパワーに関する表現であることから加算平均が可能なので、マルチルック処理を並行して行なって後方散乱係数を取得します。マルチルック処理とは、コヒーレントな電波を用いて対象に電波を送信して対象物を計測する場合、原理的に生じる画素値のノイズ(スペックル)を低減させるために行われるものです。最終的に、散乱波のパワーに関するSARデータは、対象の単位面積当たりの散乱波のパワーである後方散乱係数によって記述されます。
この過程で、SAR画像の前処理工程は完了するので、そこからシグネチャを抽出して対象の散乱の振る舞い方を解析したり、様々な画像分類アルゴリズムを使用して土地の被覆分類などを行うことが可能になります！(つまり、いろんな解析ができます笑)
実装アルゴリズム
この解析は、画像をクロップして切り出さずに解析を行なった場合(つまり、TellusのAPIからHH,HV,VV成分をそのまま取得したデータで解析を行なった場合、かなりの解析時間がかかります(ミュラー行列を算出するだけで3時間、そこからマルチルック処理が90分くらいかかるので、全体で一つのSAR撮影画像を前処理するのに5時間もかかります笑)。なので、この処理は画像を切り抜いた状態で解析を行うか、サーバ上でジョブを投げて管理しないとTellusの開発環境上では実践しにくいので注意してください汗
また、globで取り出しているのは.npyとしてバイナリデータを処理したものを使っています。SARのバイナリデータからnumpyへの変換方法は前回の記事を参考にしてください。また、今回のマルチルック処理は、アジマス方向に4ルックのマルチルック処理を施しました(つまり、縦方向に1/4に圧縮されています)。
import requests from PIL import Image import urllib.request import matplotlib.pyplot as plt import numpy as np import pandas as pd import codecs import struct from io import BytesIO from PIL import ImageOps import re import time import os import glob SAR_list = glob.glob('~/dataset/*.npy') i=0 res = [] for data in SAR_list: SAR = np.load(data) print(SAR.shape[0],SAR.shape[1],SAR.dtype) res.append(SAR) i=i+1 res = np.array(res) def multilook(im): original_width = im.shape[1] original_height = im.shape[0] width = original_width height = original_height // 4 resized_image = np.zeros(shape=(height, width), dtype=np.float64) scale = 4 for i in range(height): for j in range(width): temp = 0 for x in range(scale): for y in range(scale): temp += im[i*scale + x, j] resized_image[i, j] = temp/(scale*scale) return resized_image def mueller(HH,HV,VV,x,y): m11 = np.zeros((x,y)) m12 = np.zeros((x,y)) m13 = np.zeros((x,y)) m14 = np.zeros((x,y)) m21 = np.zeros((x,y)) m22 = np.zeros((x,y)) m23 = np.zeros((x,y)) m24 = np.zeros((x,y)) m31 = np.zeros((x,y)) m32 = np.zeros((x,y)) m33 = np.zeros((x,y)) m34 = np.zeros((x,y)) m41 = np.zeros((x,y)) m42 = np.zeros((x,y)) m43 = np.zeros((x,y)) m44 = np.zeros((x,y)) for i in range(x): for j in range(y): a = HH[i,j] b = a.conjugate() c = HV[i,j] d = c.conjugate() e = VV[i,j] f = e.conjugate() m11[i,j] = (0.25 * (a*b + 2*c*d + e*f)).real m12[i,j] = (0.25 * (a*b - e*f)).real m13[i,j] = (0.5 * ((a*d).real + (c*f).real)).real m14[i,j] = (0.5 * ((c*b).imag + (e*d).imag)).real m21[i,j] = (0.25 * (a*b - e*f)).real m22[i,j] = (0.25 * (a*b - 2*c*d + e*f)).real m23[i,j] = (0.5 * ((a*d).real - (c*f).real)).real m24[i,j] = (0.5 * ((c*b).imag + (c*f).imag)).real m31[i,j] = (0.5 * ((a*d).real + (c*f).real)).real m32[i,j] = (0.5 * ((a*d).real - (c*f).real)).real m33[i,j] = (0.5 * ((a*f).real + c*d)).real m34[i,j] = (0.5 * (e*b).imag).real m41[i,j] = (0.5 * ((c*b).imag + (e*d).imag)).real m42[i,j] = (0.5 * ((c*b).imag + (c*f).imag)).real m43[i,j] = (0.5 * (e*b).imag) m44[i,j] = (0.5 * ((c*d) - (a*f).real)).real m11 = multilook(m11) m12 = multilook(m12) m13 = multilook(m13) m14 = multilook(m14) m21 = multilook(m21) m22 = multilook(m22) m23 = multilook(m23) m24 = multilook(m24) m31 = multilook(m31) m32 = multilook(m32) m33 = multilook(m33) m34 = multilook(m34) m41 = multilook(m41) m42 = multilook(m42) m43 = multilook(m43) m44 = multilook(m44) Shh_conj_Shh = 2 * m12 + 2 * m11 - m33 - m44 Shv_conj_Shv = m33 + m44 Svv_conj_Svv = 2 * m11 - 2 * m12 - m33 - m44 Re_Svv_conj_Shh = m33 - m44 Im_Svv_conj_Shh = 2 * m34 Re_Shv_conj_Shh = m13 + m23 Im_Shv_conj_Shh = m14+ m24 Re_Shh = np.sqrt(Shh_conj_Shh ) abs_Shv = np.sqrt(Shv_conj_Shv) abs_Svv = np.sqrt(Svv_conj_Svv) lines = Re_Shh.shape[0] pixels = Re_Shh.shape[1] Shh = np.zeros( (lines,pixels), dtype=np.complex64) Im_Shh = np.zeros( (lines,pixels), dtype=np.complex64) ############### Shh.real = Re_Shh Shh.imag = Im_Shh ############### abs_Shv_conj_Shh = np.sqrt(Re_Shv_conj_Shh * Re_Shv_conj_Shh + Im_Shv_conj_Shh * Im_Shv_conj_Shh) Re_Shv = abs_Shv * Re_Shv_conj_Shh / abs_Shv_conj_Shh Im_Shv = abs_Shv * Im_Shv_conj_Shh / abs_Shv_conj_Shh Shv = np.zeros( (lines,pixels), dtype=np.complex64) ############### Shv.real = Re_Shv Shv.imag = Im_Shv ############### abs_Svv_conj_Shh = np.sqrt(Re_Svv_conj_Shh * Re_Svv_conj_Shh + Im_Svv_conj_Shh * Im_Svv_conj_Shh) Re_Svv = abs_Svv * Re_Svv_conj_Shh / abs_Svv_conj_Shh Im_Svv = abs_Svv * Im_Svv_conj_Shh / abs_Svv_conj_Shh Svv = np.zeros( (lines,pixels), dtype=np.complex64) ############### Svv.real = Re_Svv Svv.imag = Im_Svv ############### return Shh, Shv, Svv Shh, Shv, Svv = mueller(res[0],res[1],res[2],res[0].shape[0],res[0].shape[1]) np.save(os.path.join('/home/ssh_admin','Shh'),Shh) np.save(os.path.join('/home/ssh_admin','Shv'),Shv) np.save(os.path.join('/home/ssh_admin','Svv'),Svv) def SAR_plot(HH_data,HV_data,VV_data): data = np.stack([HH_data,HV_data,VV_data],axis=2) #data = 10*np.log10(data.real**2 + data.imag**2)-83-32 #-83-32はいらないです。PALSAR2の説明書にはキャリブレーションの仕様で導入されていたんですが、入れると変になります汗 data = 10*np.log10(data.real**2 + data.imag**2) amax = data.max() amin = data.min() scale = 255.0/(amax-amin) img = data*scale img = img.astype(np.uint8) img = Image.fromarray(img).convert('RGB') plt.figure(figsize=(10,10)) img.save('SAR_01.jpg', quality=95) plt.imshow(img) Shh = np.load('Shh.npy') Shv = np.load('Shv.npy') Svv = np.load('Svv.npy') SAR_plot(Shh,Shv,Svv)
結果
両者を見比べてみると、ミュラー平均化した画像の方が、画素の揺らぎ(スペックルノイズ)がかなり抑制されてめちゃくちゃ綺麗に可視化できていると思います！ このようにSAR画像は前処理の工程が非常に難しいのですが、今回実装ができたので、このSAR画像を使って
・土地の被覆分類をディープラーニングなどを応用して行う
などの様々な解析手法を適応させることが可能になるので、興味がある人はぜひこの前処理をSAR画像に適用させて様々な機械学習その他諸々の解析に応用してみてください！！
まとめ
Tellus運営さん
研究と業務委託、実習などを並列して行なってる博士後期課程の身分なので、SAR画像の解析にソースを割くのが難しい状況なのですが、今後もいろんな解析を試していく予定です。ですので、どうか何卒、開発環境の期限延長をお願い致します...!
参考資料
https://dil-opac.bosai.go.jp/publication/nied_report/PDF/63/63jitufuchi.pdf
TexのTableに番号なしの脚注を書くためのTips
LAPRASのスコアを上げたくてQiitaにネットワーク解析を投稿した結果、暫くネタがなかった...。
現在は研究結果をジャーナルとして纏めているんですが、テーブルの下に番号のない脚注の書き方にちょっと苦戦したので、Tipsとして掲載することにしました。
実現するための方法はとても簡単で、以下をプリアンブルに記述するだけです。
¥usepackage{lipsum}
¥newcommand¥blfootnote[1]{%
¥begingroup
¥renewcommand¥thefootnote{}¥footnote{#1}%
¥addtocounter{footnote}{-1}%
¥endgroup
}
実際に本文で使用する際には、以下のように¥end{table}の下に加えれば番号のない脚注を表示させることが可能です。
¥begin{table}[b]
ここにテーブルの内容を記述します。
¥end{table}
% これで番号がない脚注を書くことができます。
¥blfootnote{これでテーブルの脚注に番号が入らずに生成できる}
まとめ
今回は、Texの番号なし脚注番号を書くための方法を記載しました。そろそろちゃんとした解析を載せたい...。
参考資料
https://user.ecc.u-tokyo.ac.jp/users/user-15826/wiki/?TeX/Manual/footnote tex.stackexchange.com
PALSAR2の全偏波成分を使ってRGB形式で散乱特性を可視化する
今年は、本当に暖冬ですね。普段は関西の田舎に住んでるので、寒暖差が激しすぎて耐えられません...が、今年は大丈夫そうです笑。
前回は、Tellusと呼ばれる衛星画像のプラットフォームからAPIを使ってPALSAR2のSARデータを取り出し、取得したイメージファイル(CEOS format)のバイナリ処理をしてHH偏波の生データを取り出してグレースケールで可視化するところまでやりました。
今回は、散乱特性を可視化するために、HH偏波だけではなく全ての偏波成分(HH, HV, VV)を使って、散乱特性をRGB形式で可視化したいと思います。
RGB形式で複数の偏波成分を組み合わせて、散乱特性を可視化すると何が分かるの？
元々、偏波は「時間とともに変動する電界ベクトル先端の軌跡」を意味しています。例えば、電界が直線上に振動して伝播する偏波を直線偏波、円を描く時には円偏波などと呼びます。基本的に地表や建物など、衛星から地表へ電波が送信された時に跳ね返る対象物は様々なものがあり複雑な形状であるものが多いため、ある偏波で対象物に向かってマイクロ波を照射した場合、多くが楕円偏波としてデータが取得されます。この偏波を表現するためには少なくとも独立した二つの次元の情報が必要で、偏波レーダの場合は直線偏波である水平偏波(H)と垂直偏波(V)の二種類を組み合わせて散乱特性を表現しています。
水平偏波は、地面に水平に伝播する偏波を示していて、垂直偏波は地面に垂直に伝播する偏波で水平偏波とは直交しています。偏波レーダでは、水平偏波Hを送信して、受信された反射波の水平偏波Hと垂直偏波Vを同時に取得し、それぞれと定義することが多いです。同様に、垂直偏波を送信した場合においても、受信された反射波の水平偏波と垂直偏波を取得して
と表現されます。これら4つの成分を行列にしたものは散乱行列と呼ばれており、この散乱行列を元にして偏波情報をセンシングする技術のことを、専門分野では「レーダポラリメトリ」と呼ばれています。
前回は、マイクロ波の波長(バンド)によって見え方に違いが生じると説明しましたが、偏波成分によっても違いがあります。例えば、HH偏波は二回散乱のような建物などに対して後方散乱強度が強く、HV偏波は森林などに対して後方散乱強度が強いなどの特徴があり、これらを組み合わせてRGBで可視化する事によって、森林か建物かどうかの違いが容易に判別できます。
では、実際にSAR画像をRGB形式で可視化して見ましょう。今回は、R:HH偏波、G:HV偏波、B:VV偏波として可視化します。使用したSAR画像は山の森林地域で季節は9月の終わり頃なので、山が緑っぽく見えれば成功です。
SAR画像のイメージファイルの取得と可視化までのソースコード
まずは、必要なモジュールを読み込んで、取り出したSAR画像のイメージファイル名をHH,HV,VVの順番でソートして取得します。
import requests from PIL import Image import urllib.request import matplotlib.pyplot as plt %matplotlib inline import numpy as np import pandas as pd import codecs import struct #from skimage import io from io import BytesIO from PIL import ImageOps import re import time ## APIのドメイン URL_DOMAIN = "file.tellusxdp.com" # デフォルトのプロダクトIDでも良い場合はこちらのオプションを選択する SEARCH_DOMAIN = "api/v1/origin/search/palsar2-l11" PUBLISH_DOMAIN = "api/v1/origin/publish/palsar2-l11" # フルポラりメトリなど、使用したい衛星画像にこだわりがある場合は以下でプロダクションIDを設定する必要性がある。 #SEARCH_DOMAIN = "api/v1/palsar2/get_scene/UBSL1.1R_D" ## APIトークンをセットしてください BEARER_TOKEN = "4feb92f9-b3dc-4055-9446-a72ab6a908bf" #area = {"min_lat":35.487,"min_lon":139.394,"max_lat":35.599,"max_lon":139.736,"censor":'VV'} area = {"left_bottom_lon":137,"right_top_lon":140,"left_bottom_lat":35,"right_top_lat":37,"polarisations":"HH+HV+VV+VH"} headers = {'Authorization': 'Bearer %s' % BEARER_TOKEN} if BEARER_TOKEN == "": print("APIトークンがセットされていません") def get_data(domain, parameters, location): res = requests.get("https://{}/{}".format (domain, parameters), params=location, headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res.json() def get_URL(domain, parameters, dataset_id): res = requests.get("https://{}/{}/{}".format (domain, parameters,dataset_id), headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res.json() def get_image(url): res = requests.get(url, headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res res = get_data(URL_DOMAIN, SEARCH_DOMAIN, area) ids = res['items'][2]['dataset_id'] # 今回は、検索した三番目にヒットしたSAR画像のプロダクトを使用しました。 URL = get_URL(URL_DOMAIN, PUBLISH_DOMAIN,ids) SAR = sorted({(k['file_name'],k['url']) for k in URL['files'] if re.search('IMG-(HH|HV|VV)', str(k))}) print(list(SAR)[0][0],list(SAR)[0][1]) print(list(SAR)[1][0],list(SAR)[1][1]) print(list(SAR)[2][0],list(SAR)[2][1])
出力結果
IMG-HH-ALOS2127230710-160930-HBQR1.1__A HH偏波のURL IMG-HV-ALOS2127230710-160930-HBQR1.1__A HV偏波のURL IMG-VV-ALOS2127230710-160930-HBQR1.1__A VV偏波のURL
次に、前回使用したSAR画像の抽出方法を関数として定義しました。
def SAR_data(data): b_data = get_image(data) lines = int(b_data.content[236:244]) pixels = int(b_data.content[248:256]) real = np.zeros((lines,pixels),np.float32) imag = np.zeros((lines,pixels),np.float32) result = np.zeros( (lines,pixels), dtype=np.complex64) raw_data = b_data.content[720:] for j in range(lines): for i in range(pixels): ret = struct.unpack('>ff', raw_data[72864*(j)+544 + 8*(i):72864*(j)+544 + 8*(i)+8]) real[j,i] =ret[0] imag[j,i] =ret[1] result.real = real result.imag = imag return result, lines, pixels
実際にSAR画像のHH偏波、HV偏波、VV偏波の生データを取得します。三つのSAR画像を取得するのに、大体15分くらいかかりました(それなりにかかりました笑)
t1 = time.time()
HH_data, lines, pixels = SAR_data(list(SAR)[0][1])
HV_data, _, _ = SAR_data(list(SAR)[1][1])
VV_data, _, _ = SAR_data(list(SAR)[2][1])
# 処理後の時刻
t2 = time.time()
# 経過時間を表示
elapsed_time = t2-t1
print(f"経過時間：{elapsed_time}")
では、取得したSAR画像のそれぞれの偏波成分をRGBで可視化する関数を作って、結果を表示して見ます。 RGB形式にする場合は8or16bitに標準化する必要があるので、スケール変数を設定して各SAR画像に適用しています。
def SAR_plot(HH,HV,VV):
data = np.stack([HH_data,HV_data,VV_data],axis=2)
data = 10*np.log10(data.real**2 + data.imag**2)-83-32
amax = data.max()
amin = data.min()
scale = 255.0/(amax-amin)
img = data*scale
img = img.astype(np.uint8)
img = Image.fromarray(img).convert('RGB')
plt.figure(figsize=(10,10))
img.save('SAR_01.jpg', quality=95)
plt.imshow(img)
SAR_plot(HH_data,HV_data,VV_data)
可視化してみた結果
上図の可視化結果は全体の一部を切り取ったものですが、よくみると森林は緑色で表示されてること確認できると思いますので、正常にSAR画像の散乱特性を可視化できたのではないでしょうか(はてなブログだと形式上無理でしたが、実際のTIFFで拡大すると、細かい地域特性が見れました)。
まとめ
今回は、SAR画像の複数の偏波成分を組み合わせることで、散乱特性を反映した画像を可視化することができました。SARの散乱特性の大きな性質として、山間部のような地域では台風などによる洪水や季節変化による森林地形の変化などで散乱特性が季節や災害ごとに大きく変動しますが、都市部の環境などでは散乱特性があまり変わらないという特徴があります。これらの特性を生かして様々な時系列のSAR画像を組み合わせて様々な都市環境の特徴量を取り出すことも可能です。
次は、SARのスペックルを取り除く手法をいくつか紹介できればと思っております。
TellusのAPIからPALSAR2の生データを取得する
まず始めに
前回は、SAR(Synthetic Aperture Radar: 合成開口レーダ)に関する基礎的な知識だったり特徴などをお話しさせて頂きました。
今回は、Tellusという衛星データプラットフォームで公開されているPALSAR2というSARの画像の生データを抽出して可視化するために、
・TellusのAPIを用いてSAR画像のイメージファイルなどを取得
・イメージファイルのバイナリ情報を処理をして、位相情報が入った複素画像を取得
・実際に強度画像を表示してSARデータを取り扱ってみる。
というアプローチを行います。本当はGDALなどの便利なライブラリがpythonにはあるのですが、衛星ごとに画像のフォーマットが異なるため、既存ライブラリでは対応できない場合があります。私も最初はGDALを使ってみたのですが、公開されているPALSAR2のSAR画像を読み込むことができませんでした。なので、今回はJAXAが公開している説明書に書かれている情報を元にして必要なパラメータなどの情報を取り出し、各ピクセルごとの複素数データを抽出する必要性があります。
今回Tellusで取り扱うSAR画像について
今回取り出すのは、PALSAR2のSAR画像のプロダクトL1.1を使います(というよりも、現時点ではTellusではプロダクトL2.1が公開されていないので、GeoTIFF形式の強度画像を扱うことができません)。L1.1のプロダクトは、前回の記事で申し上げたようにCEOSという特殊なフォーマットを用いており、衛星ごとに仕様が異なります。なので、プロダクトのデータ内容を確認してバイナリデータの処理を行なってイメージデータを取り出す必要性があります。
ちなみに、最終的に可視化した結果はこちらになります。
使用した開発環境
今回は、以下の環境でJupyter Notebook上からSAR画像の生データを取り出して可視化します。TellusのサイトでJupyter環境を申請すれば誰でもAPIを通してあらゆる衛星データを扱えますので、もし興味がある方は、Tellusで開発環境を申請してみてください。ちなみに、自分は機械学習などの解析も考慮に入れているので、以下の高火力GPUコンピューティングの環境を申請して使用しております。
実際の開発環境 GPUカード：NVIDIA Tesla V100 for PCI-Express (32GB) ×4 CPU：Xeon E5-2623 v3 4Core×2 (8C/16T 3.0GHz Max3.5GHz) メモリ：128GB ディスク（RAID1）：SSD 480GB 2台/1組 Ubuntu: 18.04
高火力GPUコンピューティング環境については、実際にPytorchを導入してJupyterでリモートで操作できる環境を実装するまでの備忘録を過去の記事で公開しているので、興味がある方はこちらも参考にして環境のセットアップなどをしてみてください。
実際に取り出してみる
TellusのAPIを用いてSAR画像のイメージファイルなどを取得
SAR画像を取得するために、まずTellusのAPIを通してSAR画像のプロダクトのURL情報を抽出して、そこから取得したURLにアクセスして画像を取り出す形式になっています。TellusのAPIのリファレンスは以下で説明されているので、詳しくはこちらを参考にしてください。ソースコード上にも書いてますが、APIを利用するにはトークンが必要です。なので、Tellusにログインして、マイページにあるAPIアクセス設定で発行したAPIトークンを変数として入れる必要があります。また、APIのデフォルト設定ではHH偏波(水平偏波を送信して同じ偏波成分で受信されたもの)しか取得できないらしいです。なので、フルポラリメトリで取得したい場合には、TellusのAPIリファレンスには残念ながら載っていないのですが、"polarisations":"HH+HV+VV+VH"をパラメータとして入力する必要性があるので注意してください。
まずは、必要なライブラリとパラメータを書き込んでいきます。
#------------------------------------------------# # 必要なライブラリを読み込む import requests from PIL import Image import urllib.request import matplotlib.pyplot as plt %matplotlib inline import numpy as np import pandas as pd from io import BytesIO import codecs import struct #------------------------------------------------# #------------------------------------------------# ## APIのドメイン URL_DOMAIN = "file.tellusxdp.com" SEARCH_DOMAIN = "api/v1/origin/search/palsar2-l11" # PALSAR2の標準処理データ(L1.1)を検索する PUBLISH_DOMAIN = "api/v1/origin/publish/palsar2-l11" # データを取得するためのURLを手に入れる ## APIトークンをセットしてください BEARER_TOKEN = "Tellusにログインして、マイページにあるAPIアクセス設定で発行したAPIトークンを入れる" # 取得したい地域の緯度経度の範囲情報を入力。詳しい詳細はTellusのAPIリファレンスのページを参照してください # フルポラリメトリで衛星画像を取得したい場合は、"polarisations":"HH+HV+VV+VH"を追加します area = {"left_bottom_lon":130,"right_top_lon":140,"left_bottom_lat":35,"right_top_lat":40,"polarisations":"HH+HV+VV+VH"} headers = {'Authorization': 'Bearer %s' % BEARER_TOKEN} if BEARER_TOKEN == "": print("APIトークンがセットされていません") #------------------------------------------------#
次に、PALSAR2のデータを検索する関数、検索して出てきたプロダクトIDを元にしてSARデータを取得するためのURLを発行する関数、そのURLを元にSARデータを持ってくる関数という形で、合計で三つの関数を定義します。
# PALSAR2のデータを検索する関数(プロダクトIDを取得する) def get_data(domain, parameters, location): res = requests.get("https://{}/{}".format (domain, parameters), params=location, headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res.json() # 取得したプロダクトIDにアクセスして、SARデータを取得するためのURLを手に入れるための関数 def get_URL(domain, parameters, dataset_id): res = requests.get("https://{}/{}/{}".format (domain, parameters,dataset_id), headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res.json() # URLを元に実際のSAR画像(イメージファイル)やそのサマリーなどを取得する関数 def get_image(url): res = requests.get(url, headers=headers) if res.status_code is not 200: raise ValueError('status error({}).'.format(res.status_code)) return res
では、まずは取得したSAR画像のプロダクトIDを取得したいと思います。検索結果はデフォルトで100件出てくるのですが、今回はお試しですし多すぎるので10件だけ表示しています。パラメータとして入力した緯度経度の範囲内のプロダクトを検索してくれているみたいです。
res = get_data(URL_DOMAIN, SEARCH_DOMAIN, area) for i in range(10): print(res['items'][i]['dataset_id'])
実行結果: ALOS2127230690-160930 ALOS2127230700-160930 ALOS2127230710-160930 ALOS2127230720-160930 ALOS2127230730-160930 ALOS2127380690-161001 ALOS2127380700-161001 ALOS2127380710-161001 ALOS2127970690-161005 ALOS2127970700-161005
このようにして、実際にPALSAR2のプロダクトIDを取得することができました。
次に、対象のプロダクトIDを入れてSAR画像を取得するためのURLを手に入れます。今回は、最初に出てきたALOS2127230690-160930のプロダクトIDを入力しました。
ids = res['items'][0]['dataset_id'] SAR = get_URL(URL_DOMAIN, PUBLISH_DOMAIN,ids)
取得したURLをみると、以下のような形式になっています。
{'dataset_id': 'ALOS2127230690-160930',
'expires_at': '2020-01-09T01:20:09.452570+00:00',
'files': [{'file_name': 'BRS-HH-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 82228,
'url': 'HH偏波の縮小画像のjpegファイルを取得できるURL情報'},
{'file_name': 'BRS-VV-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 79224,
'url': 'VV偏波の縮小画像のjpegファイルを取得できるURL情報'},
{'file_name': 'IMG-VV-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': 'VV偏波のイメージファイルのURL情報'},
{'file_name': 'IMG-HH-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': ’HH偏波のイメージファイルのURL情報'},
{'file_name': 'VOL-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 2880,
'url': 'ボリュームディレクトリという、ファイルの管理情報などが格納されているファイルを取得できるURL情報'},
{'file_name': 'IMG-VH-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': 'VH偏波のイメージファイルのURL情報'},
{'file_name': 'TRL-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 568980,
'url': 'トレイラと呼ばれる、低分解能画像データが格納されているファイルを取得するためのURL情報'},
{'file_name': 'summary.txt',
'file_size': 1435,
'url': 'データのサマリ情報を取得できるURL'},
{'file_name': 'IMG-HV-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1697149008,
'url': 'HV偏波のイメージファイルのURL情報'},
{'file_name': 'LED-ALOS2127230690-160930-HBQR1.1__A',
'file_size': 1609432,
'url': 'SARリーダと呼ばれる、地図投影や姿勢、データ品質やラジオメトリックデータが格納されているファイルを取得するためのURL情報'},
{'file_name': 'BRS-HV-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 75795,
'url': 'HV偏波の縮小画像のjpegファイルを取得できるURL情報'},
{'file_name': 'BRS-VH-ALOS2127230690-160930-HBQR1.1__A.jpg',
'file_size': 75672,
'url': 'VH偏波の縮小画像のjpegファイルを取得できるURL情報'}],
'project': 'palsar2-l11'}
今回は、SARの生データ(各ピクセルに複素数データが格納されている形式のデータ)を取り出して可視化するのが目的なので、'IMG-HH-ALOS2127230690-160930-HBQR1.1__A'(つまり、今回のプログラムでいうSAR['files'][3]['url'])のURLにアクセスして、HH偏波のイメージファイルを取得します。イメージファイルはかなりデータの容量が大きいので、処理が終了するまでに少しだけ時間がかかります。
thumbs = get_image(SAR['files'][3]['url']) data = thumbs
イメージファイルのバイナリ情報を処理をして、振幅と位相情報が入った複素画像を取得する
では、TellusのAPIを通してPALSAR2のHH偏波のイメージファイルを取得できたので、次はイメージファイルのバイナリ情報を処理してデータを抽出します。これはCEOSフォーマットになっているので、実際には以下のプロダクトの詳細に目を通して細かいパラメータ情報やその位置を把握する必要があります。
https://www.eorc.jaxa.jp/ALOS-2/doc/fdata/PALSAR-2_xx_Format_CEOS_J.pdf
ここでは簡潔にいうと、取得したSARイメージファイルの最初の720バイト分はファイルディスクリプタといい、ファイル情報の詳細が載っています。そのあとに処理済データレコードがあり、この中にSAR画像の生データの情報が入っていますが、以下の図ように、画像のラインごとにSARデータの前にプレフィックスデータが544バイト分入っているので、画像のラインごとにこのプレフィックスの情報を取り除いてからデータを取り出す必要があります。(PALSAR2のCEOSフォーマットの説明書で使用されていた図を抜粋しています)
また、SARデータの後に付随しているダミーデータや不定データのバイトはデフォルトで0となっているのはいいのですが、SARデータ自体はピクセル数分だけ情報がチャンクされている形になっています。さらに、1ピクセルあたり8バイトの情報のうち、前半4バイトには2の補数形式で実数情報(振幅情報)、後半4バイトには虚数成分(位相情報)が格納されています。どちらもfloat32の形式で取り出せますが、注意するところはビッグエンディアンで取り出す必要があるところです。自分が行なったjupyter上だと、デフォルトではリトルエンディアンで取り出されるようなので、今回はstructモジュールのunpackを使ってデータを処理しているのですが、>ffの形で8バイトのデータ情報を取り出すようにしてください。
それでは、まずはファイルディスクリプタにある基本的な情報を取り出して見ましょう
print(data.content[48:64]) print(data.content[16:27]) print(data.content[181:186]) # データセットサマリレコードの数 print(data.content[236:244]) # ライン数 print(data.content[248:256]) # ピクセル数 print(data.content[186:192]) # データセットサマリレコード長 SARデータレコード長 print(data.content[296:304])
出力結果: b'AL2 SARBIMOP ' b'CEOS-SAR ' b'23292' b' 23292' b' 9040' b' 72864' b' 13 4PB'
以下のようなバイト形式で取り出せました。ここで重要なのはライン数、ピクセル数、SARデータのレコード長です。実際にPALSAR2のSARデータを取り出す時には、これらの情報が非常に重要になります。
#ライン数、ピクセル数、SARデータレコード長を変数として定義する lines = int(data.content[236:244]) pixels = int(data.content[248:256]) Length = int(data.content[186:192])
次に、必要ないファイルディスクリプタの720バイト分の情報を取り除きます。これで、raw_dataに格納されているのはラインごとにプレフィックスとピクセル数分のSARイメージ情報が入ったデータのみとなります。
raw _data= data.content[720:]
試しに、最初のラインに入っているプレフィックスデータの情報を一部見て見ましょう。プロダクトの説明書通りであれば、左と右詰めのデータはないはずなので、ダミーデータなどはなく純粋にプレフィックスデータとSARデータのみで構成されていることが確認できると思います。
print(int(codecs.encode(raw_data[20:24], 'hex') ,16)) # 左詰め*(デフォルトは0) print(int(codecs.encode(raw_data[24:28], 'hex') ,16)) # ピクセル数 print(int(codecs.encode(raw_data[28:32], 'hex') ,16)) # 右詰め*(デフォルトは0) print(int(codecs.encode(raw_data[48:50], 'hex') ,16)) # フルポラリメトリのプロダクトなら4が出てくる
出力結果: 0 9040 0 4
では、実際にSARの生データをバイナリ情報を処理して抽出して見ます。
real = np.zeros((lines,pixels),np.float32) imag = np.zeros((lines,pixels),np.float32) import time t1 = time.time() for j in range(lines): for i in range(pixels): ret = struct.unpack('>ff', raw_data[Length*(j)+544 + 8*(i):Length*(j)+544 + 8*(i)+8]) real[j,i] =ret[0] imag[j,i] =ret[1] # 処理後の時刻 t2 = time.time() # 経過時間を表示 elapsed_time = t2-t1 print(f"経過時間：{elapsed_time}") # 自分の環境では、大体三分ちょいで処理が終わりました。
これで、ようやく各ピクセルに複素数の情報が格納されたSAR画像を取得することができました。
実際に強度画像を表示してSARデータを取り扱ってみる。
では、最後に複素数の情報を強度画像に変換して、グレースケールで可視化します。
from PIL import ImageOps result4 = np.zeros( (lines,pixels), dtype=np.complex64) result4.real = real result4.imag = imag resres = 10*np.log10(np.abs(result4)) print(resres.max(),resres.min()) img = Image.fromarray(resres).convert('L') plt.figure(figsize=(16,16)) plt.imshow(img, cmap='gray', interpolation='nearest')
可視化した結果が、前回の記事で例としてお見せしたSAR画像になっていて、以下の画像になります。
これで、ようやくPALSAR2の生データの取り出しと可視化に成功しました。これ以降は抽出したこれらの画像に対してスペックルのキャリブレーションを行なったり、他の偏波成分(HVやVV偏波)の複素数データを取り出して、HH,HV,VVの順番でRGB形式にしてその地域の散乱特性情報を可視化したりなど、様々なデータの加工や解析ができます。
まとめ
今回は、TellusのAPIを通じてPALSAR2の複素数データ(生データ)を取得して可視化するまでの一連の手順を紹介させていただきました。今後は、補正処理や実際のSAR解析の一端をご紹介できればと思っておりますので、よろしくお願い致します。
もしもpythonではなく、他の言語でSAR画像を研究などの都合で取り出したい場合は、昨日探していて見つけた記事ですが以下のものがありましたので、こちらもぜひ参考にして見てください(最近の記事として投稿されていたので驚きました笑)
(SAR画像解析をする前に)SARってなんぞ？
まず始めに
久しぶりの更新になってしまいました。今回の風邪は長引くのが特徴らしく、珍しく12月中旬に二週間近く風邪を引いて、また年始に五日間くらい風邪を引くハプニングにあいました。皆さん風邪には気をつけてください。
現在、趣味の範囲ですがTellusのGPU高火力サーバの環境を使って衛星データなどの解析をしております。
私はSAR画像に興味があるので、上記の環境を使って取得したSAR画像を解析して何らかの発見ができればと思っているのですが、SAR画像ってあまり使われるイメージがないので、Google Earthでよくみるような画像とどのような違いがあるのか分からない方も多いのではないでしょうか。
そこで、今回はSAR画像の解析などをブログで公開する前に、一旦ここでSAR画像についてお話できればと思います。
SARって何？
まず、SARは合成開口レーダ(Synthetic Aperture Radar)の略称で、環境変化などを地球規模で観測することができるリモートセンシング技術の一つです。最大の特徴は、何と言っても昼夜天候を問わず観測することが可能という点です。通常、光学センサで観測するライダは、雲がある場合は貫通することができず観測することができないですが、SARはマイクロ波を使って観測しており、マイクロ波は雲も貫通するので観測が可能です。つまり、災害時に雲で覆われていたり、火事で煙が覆っていたとしても観測できるなどの利点があります。また、乾いた土であるなどの条件はありますが、土の中を潜って観測することができるので、土壌や建物内部の情報も観測することができます。
さらに、SAR画像の観測に用いられているマイクロ波は、波長によってPバンド、Lバンド、Cバンド、Xバンドに分類されます。実は、これらの波長によって地表へ反射する特性が大きく異なります。SARは地表に向けて電波を送信し、対象物から反射された電波を受信していますが、この対象物が例えば表面が荒いかどうかなどの特性によって、短い波長であれば表面の凸凹に影響され、波長が長いと影響されない特徴があります。一番分かりやすい例は対象物が森林の場合で、波長が短いと樹冠から反射しますが、波長が長くなればなるほど木の幹の中で多重散乱したり、木の根元から反射したりします。実際、自分自身もこの反射特性の違いを利用して、森林に覆われた川などを可視化した経験があります。
主に使用されているバンドは国ごとに異なっていて、例えばカナダなどのように森林が多い地域では、森林の監視などの用途に用いることからCバンドが利用される場合が多いようです。日本では、Lバンドが主に使用されています。代表的なのはPALSAR2で、衛星データのオープンデータが公開されているプラットフォームであるTellusにもPALSAR2のデータが公開されています。
日本のSAR技術は実はバカにならなくて、一般的にSARというとPOLSAR(Polarimetric SAR)を意味していて、偏波情報を利用して観測を行なっているのですが、日本のPALSAR2の前身でALOS「だいち」に搭載されたPALSARは、世界で初めてフルポラリメトリック観測(4つの偏波情報全てを利用した観測手法)を可能にしています。
一方で、SAR画像のデメリットとしては、一般的なデジカメなどで利用されている光学センサなどとは異なり、複数のバンドを組み合わせて馴染みのあるRGB画像にすることはできず、基本的に白黒画像です。もちろん、フルポラリメトリなどと言って、全ての偏波情報を組み合わせてカラー画像にするやり方はあるのですが、各色の情報は「色」ではなく「散乱特性」という解釈をする必要があるので、あまりSAR画像に馴染みのないひとは解釈に苦労します。また、このあと実際の画像でお見せしますが、画像に白い斑点のようなものができてしまい、これらのノイズを取り除く必要もあります。さらに、山岳地帯の倒れこみが激しかったり山並みが重なって平らに見えたりなどの欠点があります。
実際のSAR画像
SAR画像って言ってもピンとこないと思うので、先日取り出したSAR画像をお見せしたいと思います(可視化方法については後日ブログで発進させて頂きます)。
どこらへんか調べてないのですが、おそらく湖とその周辺の山だと思います笑 この画像、よくみると白いモヤモヤみたいなのが見えるかなと思います。これが先ほど説明したSAR画像特有の白い斑点で、「スペックル」と言います。具体的には、画像内における各ピクセルの反射波(明るさ)に関連した統計的な変動あるいは不確定要素のことを意味します。SAR画像で解析を行う前には、ほぼ必ずと言っていいほどこのスペックルの前処理をします。よくやるのは、分解能が犠牲になりますがマルチルック(隣接するピクセルを平均化してラジオメトリック分解能を小さくする)や乗法モデルのスペックルを仮定して局所統計量を利用したスペックルフィルタリングによるアプローチがあります。他にも、SAR画像の前処理方法としては、レーダ装置が持つ偏波特性の誤差(アンテナの不完全性や宇宙環境・打ち上げの時の外的な要因)を抑えるポラリメトリック補正や散乱値(散乱特性:送信パルスと受信パルスの電力比率の変動を最小限にする処理)の変動を最小限にするラジオメトリック補正などがあります。
SARの生データってどんな感じ？
一般的なSARの生データ(PALSAR2のL1.1プロダクト)は、スラントレンジ上(センサとターゲットを結んだ視線方向)に等間隔に配置された複素数データ(SLC: Single-Look Complex)の情報が含まれていて、各ピクセルに実数(振幅)と虚数(位相)の情報が格納されている珍しい画像です。可視化する時は強度に変換する場合が多いです。
複素数データって何に使えるの？
SARのデータの性質上、一番の特徴は複素数情報があることがです。特に、位相情報が非常に有用で、例えばこの情報を利用したものにインターフェロメトリ解析があります。インターフォロメトリにもいくつか種類がありますが、一番有名なのは差分干渉SARで、これは二つのSARデータとDEMと言われる数値標高モデルを用いて、地震などによる地盤変動などの観測ができます。また、これらの位相情報を用いた位置合わせ方法もあり、二つの画像をフーリエ変換して、位相情報に変換して最も相関係数の高いピクセルの座標から二つの画像の平行移動量、さらに対数や極座標によって拡大縮小・回転角度も算出することができます。また、SARにも様々な分解方があり、偏波特性から散乱特性ごとに寄与している成分を抽出して欲しい情報を可視化する方法もあります。確か、手法によっては橋の検出なども可能だったかと思います。
SARでよく使われるデータフォーマット
SAR画像でよく使用されるフォーマットはGeoTiffとCEOSです。GeoTiffはある程度推測ができると思うので割愛しますが、ほとんどのSAR画像で複素数情報が入っている場合に関してはCEOSフォーマットになっていると思います。
CEOSフォーマットは特殊なフォーマットなので衛星ごとに仕様が異なります(なので、衛星のプロダクトの説明書を読み込む必要があります汗)。説明は以下のサイトがわかりやすいので割愛しますが、基本はデータの詳細であるイメージディスクリプタに続いて画像情報、プレフィックスなどの情報が画像の列(ライン)ごとに並んでおり、実際の画像情報を取り出す際にはファイルの中の画像情報だけを正確に取りだす必要があります。別にpythonでもC++でもMatlabでも情報は取り出せますが、バイナリデータの処理に慣れてないと大分苦労すると思います。
https://www.opengis.co.jp/htm/gamma/htm/gamma_mail_00012.htm
まとめ
今回は、SAR画像解析をする前段階としてSARについて利点や特徴などを説明しました。ここの説明だけだと偏波などの部分は分からないと思うので、もし良かったら以下のURLを参考にしていただければと思います。次回は、いよいよTellusのAPIから取得したSARイメージファイルを可視化したいと思います。よろしくお願い致します。
GPUサーバからリモートでjupyter操作してtorchも触れるようにする
前回の記事に続いて、Tellusで申請したGPUサーバーの環境にjupyterを入れて遠隔で操作できるようにし、かつpytorchも無事に使える状態にします。環境のセットアップの際にごちゃごちゃしてしまったので、もしかしたら手順が間違っている可能性がありますが、一番つまづいた根本的な原因は分かっているので、そこさえ気をつければ環境構築は可能だと思います。
前回の記事はこちらです。
そして、再喝ですがGPUのスペックは以下の通りです。
- GPUカード：NVIDIA Tesla V100 for PCI-Express (32GB) ×4
- CPU：Xeon E5-2623 v3 4Core×2 (8C/16T 3.0GHz Max3.5GHz)
- メモリ：128GB
- ディスク（RAID1）：SSD 480GB 2台/1組
- Ubuntu 18.04
Jupyter 環境の構築
まず、jupyterlab勢なので、このコマンドをします。
conda install -c conda-forge jupyterlab
Jupyterカーネルを切り替えるパッケージの設定
次に、Jupyterでカーネル(作成した仮想環境)を切り替えるためのパッケージをインストールします。
pip install environment_kernels
そして、以下のようにしてjupyterの設定ファイルを生成します。
jupyter notebook --generate-config
おそらく、前回の記事と同じ環境で設定をしているのであれば、/home/(ログインID)/.jupyter/jupyter_notebook_config.pyの方に設定ファイルができているはずなので、その中のファイルの中に以下の文を挿入します。
c.EnvironmentKernelSpecManager.conda_env_dirs = [ 'pytorchが入っている仮想環境へのファイルパス' ]
前回の記事と同じ手順を踏んで入れば、上記に入れるファイルパスは、
/home/(ログインID)/.conda/envs/(仮想環境名)/
で問題ないかと思います。
しかしながら、このままだとjupyterを起動しても(あとで繋げ方は載せてあります)、~/anaconda3/share/jupyter/kernelsに作った仮想環境のディレクトリがないので認識されません。なので、~/.conda/envsの中に入っているはずの作成した仮想環境を写す必要があります。なので、まずは~/anaconda3.share.jupyter/kernelsの下に仮想環境をコピーしてください。
そして、コピーが終わったら、次に~/anaconda3/share/jupyter/kernels/(作成した仮想環境名)の中にkernel.jsonと呼ばれるものがあるので、それを開いてみてください。
一行目が/home/(ログインID)/anaconda3/bin/pythonとありますが、このままだとデフォルトのpythonが呼び出されてしまうため、jupyterを起動してimport torchとしてもモジュールを呼び出すことができません。なので、この行を以下のように変更すれば大丈夫です。ここが色々やっても動かなかった根本的な原因でした。ウェブにも探してもなくて、前に投稿したブログを元に思い出しました。
"/home/(ログインID)/.conda/envs/(仮想環境名)",
これでjupyterを起動しても、無事にtorchなどはインポートされると思います。
Jupyter の起動方法
次に、Jupyterの起動方法を説明します。まず、GPUサーバにログインしたあと、以下のコマンドを打ちます。
jupyter lab --no-browser --port=8080
そして、別のタブのターミナルを開いて、以下を打ちます
ssh -N -L 8080:localhost:8080 (ログインID)@(環境ホスト名)
パスワードを求められるので、Tellusを利用しているのであればマイページにあるトークン情報を入力して、そのあとブラウザでhttp://localhost:8080と検索すれば、jupyterに移動できます。最初はトークンを入力するように言われるので、GPUサーバでjupyterを起動した際に表示されたURLのtoken=から続く暗号のようなものをコピペして貼り付ければjupyterの環境をリモートで操作することができるようになります。
まとめ
今回は、TellusのGPUサーバからリモートでJupyter環境を使えるセットアップを説明しました。
