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のスペックルを取り除く手法をいくつか紹介できればと思っております。
