海老ノート

Google Earth Engine 苦闘の記録

SARの勉強3:スペックル その2

フィルタリングの効果を見てみたい

ちょっと前に公開されたGoogle Earth Engineで使えるSentinel-1の前処理法を使ってフィルタリングする。以前の続きです。 前処理法はこちら https://github.com/adugnag/gee_s1_ard

フィルタリングの方法は5種類

スペックルノイズを除去する前処理法は色々あって、今でも新しい方法が提唱されたりしているようですが、ここで紹介されているのは5種類。それぞれ1枚の画像だけをつかうのか、複数枚で時間方向にもフィルタをかけるのか?で5x10種類の処理が可能。

  • Boxcar フィルタ
  • Gamma maximum a-posterior (MAP) フィルタ
  • Lee フィルタ
  • Improved Lee Sigma speckle フィルタ
  • Refined Lee speckle フィルタ

それぞれの方法で一長一短あるようで、検索すると比較する論文とか出てきますが、まあ分析対象の場所に寄りけりということなんでしょう。ざっと見た感じ最近ではImproved Lee Sigma が多く使われてるよう。ちなみにこの前処理法の論文(Mullisa et al 2021)では自分で調べてねというようなことが書いてあります。

Users can evaluate the performance of speckle filters qualitatively by visually checking the speckle reduction in homogeneous regions and the preservation of subtle features, such as roads and point scatterers, in the filtered images or quantitatively using the equivalent number of looks (ENL) or coefficient of variation in homogeneous regions in the image.

ほかに決めること
フィルタリングの窓サイズ

GEEのee.Kernelを使うので正方形のn x nの窓ではなく中心のピクセルから直径nピクセルの円として指定(数字は奇数のみ)。とりあえず5としてみました。

時間方向にフィルタリングするとき使用する画像の数

8枚分を使うとしました。これは前後まとめて3か月分のデータくらいでやってみようという見積もり。

どうやって比べるか

引用した論文の通り、画像を見た眼でみるのと均質な場所での統計量をみます。地形補正(Terrain Flattening)済みの画像を基準として、処理したものと比較。

フィルタリングした

画像をみる
ターミナル付近

この辺。

まずは1枚を処理したやつ。R:VV、G:VH、B:VV/VHで表示。

f:id:camarao:20210921131416j:plainf:id:camarao:20210921132058j:plainf:id:camarao:20210921132104j:plain
左:フィルタなし、中:Boxcar、右:GMA

f:id:camarao:20210921132107j:plainf:id:camarao:20210921132113j:plainf:id:camarao:20210921132119j:plain
左:Lee、中:Improved Lee Sigma、右:Refined Lee

時間方向にもフィルタリングすると

f:id:camarao:20210921131416j:plainf:id:camarao:20210921132438j:plainf:id:camarao:20210921132441j:plain
左:フィルタなし、中:Boxcar、右:GMA

f:id:camarao:20210921132447j:plainf:id:camarao:20210921132454j:plainf:id:camarao:20210921132459j:plain
左:Lee、中:Improved Lee Sigma、右:Refined Lee

東側のゴルフ場

1枚を処理

f:id:camarao:20210921132719j:plainf:id:camarao:20210921132724j:plainf:id:camarao:20210921132732j:plain
左:フィルタなし、中:Boxcar、右:GMA

f:id:camarao:20210921132741j:plainf:id:camarao:20210921132747j:plainf:id:camarao:20210921132753j:plain
左:Lee、中:Improved Lee Sigma、右:Refined Lee

時間方向にフィルタ

f:id:camarao:20210921132719j:plainf:id:camarao:20210921132856j:plainf:id:camarao:20210921132902j:plain
左:フィルタなし、中:Boxcar、右:GMA

f:id:camarao:20210921132907j:plainf:id:camarao:20210921132914j:plainf:id:camarao:20210921132935j:plain
左:Lee、中:Improved Lee Sigma、右:Refined Lee

滑走路南端

ここだけ色の値を変えています(値が低くて真っ暗になるから)。バンドの割り振りは同じ。 まずは1枚。

f:id:camarao:20210921133202j:plainf:id:camarao:20210921133209j:plainf:id:camarao:20210921133215j:plain
左:フィルタなし、中:Boxcar、右:GMA

f:id:camarao:20210921133223j:plainf:id:camarao:20210921133229j:plainf:id:camarao:20210921133235j:plain
左:Lee、中:Improved Lee Sigma、右:Refined Lee

取得時期の異なる複数枚でフィルタリング。

f:id:camarao:20210921133202j:plainf:id:camarao:20210921133354j:plainf:id:camarao:20210921133402j:plain
左:フィルタなし、中:Boxcar、右:GMA

f:id:camarao:20210921133409j:plainf:id:camarao:20210921133415j:plainf:id:camarao:20210921133428j:plain
左:Lee、中:Improved Lee Sigma、右:Refined Lee

感想
  • フィルタをかけるとやっぱりぼやける。ぼやけた印象が強い順にBoxcar > Lee > Improved Lee Sigma>GMA >=Refined Lee > 処理なし。

  • フィルタをかけるとノイズは減る。ノイズの少ない印象順にBoxcar > Lee > Improved Lee Sigma > GMA >=Refined Lee > 処理なし。

  • 複数枚を使った場合、ノイズは増えてぼやけた印象は減る。細部の様子がよりはっきり見える。

  • 理由は分からないけどGMAは画像を書き出すと変に切り取られた。ちょっとずれた場所なので要注意。

平均とか分散とか

とりあえずまとめ

f:id:camarao:20210921141953p:plain
CVのグラフ
見た目の評価と傾向は同じ。

  • フィルタリングするとCV(変動係数=標準偏差/平均)は大体下がる。下がり幅(=ノイズの減り)が大きいのはLeeやImproved Lee Sigma、反対にあまり下がらないのはRefined Lee。

  • 複数枚を使うとCVは上がる。

  • 元々CVの高い滑走路_VVはフィルタリングすると処理前よりCVが上がることがある。これは平均的な値をとる=ノイズに引っ張られちゃうっていうことかな。

ヒストグラムも書いてみた。

まずは処理前。

f:id:camarao:20210921142531p:plainf:id:camarao:20210921142534p:plain
左:VV、右:VH

まずはBoxcar。処理前の分布は薄い色で表示。

f:id:camarao:20210921142724p:plainf:id:camarao:20210921142728p:plain
左:VV、右:VH

GMA

f:id:camarao:20210921142912p:plainf:id:camarao:20210921142915p:plain
左:VV、右:VH

Lee

f:id:camarao:20210921142943p:plainf:id:camarao:20210921142956p:plain
左:VV、右:VH

Improved Lee Sigma

f:id:camarao:20210921143530p:plainf:id:camarao:20210921143533p:plain
左:VV、右:VH

Refined Lee

f:id:camarao:20210921143611p:plainf:id:camarao:20210921143618p:plain
左:VV、右:VH

まあ、こんな感じでフィルタをかけると値のばらつきが減って山が尖るということなんだね。グラフだらけになるので複数枚のグラフは載せないけど、全体として複数枚の方が山は平たんになります。

いったんまとめ

  • フィルタリングするとノイズは減ってエッジも減る。

  • 程度はあるけど、ノイズだけを減らしてエッジを残すというのはできない。どっちも一緒に減る。

  • なので、自分の対象とする場所の特性に合わせて丁度いいバランスのフィルタ(とその設定)を探すのがポイント。

  • 例えば分類に使うのであれば、フィルタリングした後のヒストグラムの分布の重なり具合をみて調整するとかできそう。必要以上に強いノイズ除去はせずにとどめておくような加減が探れそう。

  • 想像するに、いろんな細かい土地被覆が複雑に分布しているところはエッジ重視。反対に、スケールの大きな、一つの土地被覆パッチが大きいところではノイズを減らすのに重きを置くといい気がする。

  • 林縁とか土地被覆の境目(特にそこに道路が通ってたりすると)、反射の割合が変わるから突然周りと値が全然違う場所ができることがある(衛星の角度による)。なので、農地や森林の代表データを取りたいというときはそういう端っこを避けて取った方が綺麗なデータがとれるはず。

  • 今回は建物を入れるのを忘れてた。

  • Google Earth Engine上でフィルタリングした画像データをAssetsにエクスポートする処理がめっちゃ時間かかる。今回のサイズで2時間程度。Googleドライブへの書き出しだと2分くらいで終わるのに…多分コンピューターも頑張らないとできない処理なんだと思う。

これだけ色々やってようやく雰囲気がつかめてきた。今回は窓サイズを変えていないけど、そっちも調整できます。

コード置き場

前回に続いて慣れないGoogle Colabを使ってPythonを書きました。少しコードを書き残しておく。相変わらずよくわかんない。

注目するポイントの周りの画像を切り出して保存する処理。画像をどう取りに行けばわからなかったので、requestsでお茶を濁したけどもっといい方法があるんだと思う:

from IPython.display import Image, display_jpeg
import requests

#見たい場所の座標指定。.buffer()の引数mで範囲指定。
poi_1 = ee.Geometry.Point(141.6837, 42.7878).buffer(500) # near the terminal 500m
poi_2 = ee.Geometry.Point(141.7278, 42.7741).buffer(500) # golf course
poi_3 = ee.Geometry.Point(141.69298, 42.76525).buffer(500) # south end of a runway

#for loop の仕込み用リスト
pois = [poi_1, poi_2]
pois_title = ["poi_1", "poi_2"]

prep = [flattened, monoBoxcar, monoGamma, monoLee, monoLeeSigma, monoRefinedLee]
prepro_title = ["flattenOnly", "monoBoxcar", "monoGamma", "monoLee", "monoLeeSigma", "monoRefinedLee"]

#メインの関数
def toRgbImage(img, roi, file_name):
  rgbImage = ee.Image.rgb(img.select('VV'),
                          img.select('VH'),
                          img.select('VV').divide(img.select('VH'))) #バンドのRGBへの割り振り

  #画像が書き出されるurlを参照
  url = rgbImage.getThumbUrl({
    'min': [0.01, 0.0032, 1.25],'max': [1.0, 0.31, 31.62], 'dimensions': [512,512], 'region': roi, 'format':'jpg'}) #poi_1,2の場合の色指定
  #  'min': [0.01, 0.0032, 1.25],'max': [0.1, 0.01, 31.62], 'dimensions': [512,512], 'region': roi, 'format':'jpg'}) #poi_3の場合の色指定
  thumb = Image(url = url) 

  #urlからGoogle Colabに画像を持ってくる
  response = requests.get(url) 
  outputImg = response.content
  with open(file_name, "wb") as aaa:
    aaa.write(outputImg)
  files.download(file_name)  #自分のPCにダウンロード

  return thumb

#ループで各データ、各地点の画像を取得する
for i in range(len(prep)):
  for j in range(len(pois)):
    toRgbImage(prep[i], pois[j], prepro_title[i]+'_'+pois_title[j]+'.jpg')

重ね合わせたヒストグラムを書く:

from scipy.interpolate import make_interp_spline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

%matplotlib inline

#フィルタ処理済みの使うデータのリスト
prep = [flattened, monoBoxcar, monoGamma, monoLee, monoLeeSigma, monoRefinedLee]
prepro_title = ["flattenOnly", "monoBoxcar", "monoGamma", "monoLee", "monoLeeSigma", "monoRefinedLee"]

#注目する土地被覆のポリゴン
lc = [aoiGrass, aoiForest, aoiPavement]
lc_title = ['GrassLand', 'Forest', 'Pavement']

#線の色
color = ['darkorange', 'gold','darkgreen', 'greenyellow', 'red', 'salmon'] 

# ヒストグラムを書く関数
def combHist(prepType, pol):
  y_df = pd.DataFrame()
  for i in range(len(lc)):
    hist = prepType.select(pol).reduceRegion(
       ee.Reducer.fixedHistogram(0, 0.55, 500),lc[i]).get(pol).getInfo() # fixedHistogram(min, max, steps)
    a1 = np.array(hist)
    x1 = a1[:, 0]                 # array of bucket edge positions
    y1 = a1[:, 1]/np.sum(a1[:, 1]) # normalized array of bucket contents  
    y_df[lc_title[i]] = y1

    hist2 = flattened.select(pol).reduceRegion(
       ee.Reducer.fixedHistogram(0, 0.55, 500),lc[i]).get(pol).getInfo() # fixedHistogram(min, max, steps)
    a2 = np.array(hist2)
    x2 = a2[:, 0]                 # array of bucket edge positions
    y2 = a2[:, 1]/np.sum(a2[:, 1]) # normalized array of bucket contents  
    y_df[lc_title[i]+'_flattened'] = y2

  y_df['x'] = x1
  return y_df

# リストにしたデータごとにループ処理でグラフを書く
for k in range(len(prep)):
  plt.figure(figsize=(10,8))
  y_df = combHist(prep[k], 'VV') 
  for j in range(6):
    model=make_interp_spline(y_df['x'], y_df.iloc[:, [j]], 3)
    xs=np.linspace(0,0.5,100)
    ys=model(xs)
    plt.plot(xs, ys, c = color[j], label = y_df.columns[j] )
    plt.title(prepro_title[k]+'_VV')
    plt.xlim(0, 0.5)
    plt.legend()
  plt.savefig(prepro_title[k]+'_VV.png')
  files.download(prepro_title[k]+'_VV.png')


for k in range(len(prep)):
  plt.figure(figsize=(10,8))
  y_df = combHist(prep[k], 'VH') 
  for j in range(6):
    model=make_interp_spline(y_df['x'], y_df.iloc[:, [j]], 3)
    xs=np.linspace(0,0.12,200)
    ys=model(xs)
    plt.plot(xs, ys, c = color[j], label = y_df.columns[j])
    plt.title(prepro_title[k]+'_VH' )
    plt.xlim(0, 0.12)
    plt.legend()
  plt.savefig(prepro_title[k]+'_VH.png')
  files.download(prepro_title[k]+'_VH.png')