SAR Handbook[Chp3-2]を写経し、難しい単語を使わず意訳してみた

はじめに

本記事シリーズではSERVIRが以下のサイトで公開している「SAR Handbook」を写経し、理解を深めることを目的としています。
https://servirglobal.net/resources/sar-handbook
もともとのサイトや資料は英語にて記載されているので、翻訳しつつコード理解をしていきます。

以下の記事の続きで、本記事ではChapter3のPart2について記載します。
https://qiita.com/oz_oz/items/81ebf76aea2d930dee49

CHAPTER3

Using SAR Data for Mapping Deforestation and Forest Degradation

Chp3の説明資料は以下にあります。
https://gis1.servirglobal.net/TrainingMaterials/SAR/Ch3-Content.pdf
Chp3では主にSARによる森林監視について触れており、Sentinel-1などのSARミッションの登場により、適切な前処理と変化検出手法を用いることで、中解像度での森林監視が実現可能となったことを説明しています。

SAR Training Workshop for Forest Applications

Chp3の問題資料は以下にあります。
https://gis1.servirglobal.net/TrainingMaterials/SAR/3B_TrainingModule.pdf
Jupyter Notebook形式の資料は以下にあります。
https://gis1.servirglobal.net/TrainingMaterials/SAR/Scripts_ch3.zip

PART 2 - SAR TIME SERIES VISUALIZATION AND ANIMATIONS

Josef Kellndorfer, Ph.D., President and Senior Scientist, Earth Big Data, LLC
Revision date: January 2018
This section introduces more sophisticated animations for time series visualization which allow us to insplect time series in more depth. Note that html animations are not exported into the pdf file, but will display interactively

このPartでは、時系列変化を理解しやすくするため、アニメーションを表示する方法を紹介しています。

まずはPart1と同様にライブラリをimportします。

# Turn on inline presentations
%matplotlib inline

```python
# Imports
import os
import time
import gdal
import pandas as pd

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches  # Needed to draw rectangles
from matplotlib import animation, rc
from IPython.display import HTML

今回はアニメーション表示を行うため、matplotlibから複数のライブラリをimportしています。

データの取得方法もPart1と同様です。

West Africa - Biomass Site

datadirectory='/Users/rmuench/Downloads/wa/BIOsS1'
datefile='S32631X398020Y1315440sS1_A_vv_0001_mtfil.dates'
imagefile='S32631X398020Y1315440sS1_A_vv_0001_mtfil.vrt'
imagefile_cross='S32631X398020Y1315440sS1_A_vh_0001_mtfil.vrt'

West Africa - Niamey Deforestation Site

# datadirectory='/dev/shm/projects/c401servir/wa/cra/'
# datefile='S32631X402380Y1491460sS1_A_vv_0001_A_mtfil.dates'
# imagefile='S32631X402380Y1491460sS1_A_vv_0001_A_mtfil.vrt'

West Africa - Dam Site

# datadirectory='/dev/shm/projects/c401servir/wa/DAMsS1/'
# datefile='S32631X232140Y1614300sS1_A_vh_0001_A_mtfil.dates'
# imagefile='S32631X232140Y1614300sS1_A_vh_0001_A_mtfil.vrt'

HKH Site

# datadirectory='C:/data/hkh/time_series/S32644X696260Y3052060sS1-EBD'
# datefile='S32644X696260Y3052060sS1_D_vv_0092_mtfil.dates'
# imagefile='S32644X696260Y3052060sS1_D_vv_0092_mtfil.vrt'
# imagefile_cross='S32644X696260Y3052060sS1_D_vh_0092_mtfil.vrt'

Prepare the Animations

日付ファイルから時系列の日付を読み込み、pandasにて日付インデックスを作成します

os.chdir(datadirectory)

# Get the date indices via pandas
dates=open(datefile).readlines()
tindex=pd.DatetimeIndex(dates)

tindex

gdal.Open()関数を使用して画像を読み込みます。
また、読み込んだ画像は、上で表示した日付の画像をbandとして持っているため、.GetRasterBand()関数を使用して、1枚の画像を取得します。
今回は(1)を指定しているため、2番目の日付の'2015-04-03'の画像を取得しています。

subsetは指定しない場合は、画像が全て表示されるように設定されています。

# Open the image and read the first raster band
img = gdal.Open(imagefile)
band = img.GetRasterBand(1)
# Set the subset
if subset==None:
subset=(0,0,img.RasterXSize,img.RasterYSize)

subset

取得した画像を表示します。
画像を見やすくするため、ピクセル値の5％〜95%に絞って表示することで、極端な値の影響を抑えて、データの中心的な部分を強調することができます。

# Plot one band and subset outline to see which subset we are interested in
raster=band.ReadAsArray()
vmin=np.percentile(raster.flatten(),5)
vmax=np.percentile(raster.flatten(),95)
fig=plt.figure(figsize=(10,10))
ax=fig.add_subplot(111)
ax.imshow(raster,cmap='gray',vmin=vmin,vmax=vmax)
# plot the subset as rectangle
_=ax.add_patch(patches.Rectangle((subset[0],subset[1]),subset[2],subset[3],
fill=False,edgecolor='red'))

コードで実行していることを簡単にリスト化すると以下です。

• 画像の配列化
  • 'band.ReadAsArray()'で画像を配列で読み込んでいます。
• 画像の表示範囲の算出
  • 'raster.flatten()'で画像データを1次元にします。
  • 'np.percentile()'で任意のパーセントの位置にある値を算出し、vmin/vmaxを設定しています。
• 画像の表示
  • 'ax.imshow(raster,cmap='gray',vmin=vmin,vmax=vmax)'でvmin/vmaxを指定して、表示する画像の範囲を絞っています。
• subset範囲の強調
  • '_=ax.add_patch()'の部分でsubsetの範囲に赤い矩形を表示して、指定した範囲を強調しています。
  • 今回はsubsetを全体に設定しているのでグラフの外枠が赤くなっています。

subset領域を切り抜いた画像データを取得します。

• rastar0は1枚の画像のデータです。
• rasterstackは複数の画像のデータです。

raster0 = band.ReadAsArray(*subset)
bandnbr=0 # Needed for updates
rasterstack=img.ReadAsArray(*subset)

アニメーションを表示します。

%%capture
import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np

fig=plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.axis('off')
vmin=np.percentile(rasterstack.flatten(),5)
vmax=np.percentile(rasterstack.flatten(),95)
im = ax.imshow(raster0,cmap='gray',vmin=vmin,vmax=vmax)
ax.set_title("{}".format(tindex[0].date()))

def animate(i):
    ax.set_title("{}".format(tindex[i].date()))
    im.set_data(rasterstack[i])

# Interval is given in milliseconds
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=rasterstack.shape[0], interval=400)

コードで実行していることを簡単にリスト化すると以下です。

• 'np.percentile()'で画像の表示範囲を絞るコードは上述の内容と同じです。
• animate()関数を定義しています。
  • i番目の画像とdateを取得して、表示する画像にsetしています。
• 'matplotlib.animation.FuncAnimation()'を使用してアニメーションデータを作成しています。

matplotlibのアニメーションの埋め込み制限を変更します。
'embed_limit' はアニメーションの埋め込みに関する制限を設定します。この制限はバイト単位で表され、アニメーションを埋め込むためのメモリ制限を示します。通常、この制限は大きなアニメーションファイルを生成する場合に適用されます。
デフォルトでは、matplotlibはアニメーションの埋め込みサイズの制限を設定しています。しかし、このコードではその制限を大きな値に設定しています。これにより、より大きなアニメーションを生成しても、制限によって途中でアニメーションが切れることがなくなります。

rc('animation',embed_limit=40971520.0) # We need to increase the
# limit to show the entire animation

アニメーションを表示します。

from IPython.display import HTML
HTML(ani.to_jshtml())

Plot the global means of the Time Series for the Subset

1. Conversion to power
2. Compute means
3. Convert to dB
4. Plot time series of means

subsetの時系列平均値をプロットします。

SARの画像データは、元々の振幅値が線形にスケールされた形で保存されています。これらのデジタル数値（DN）は、適切な後方散乱値(地表や物体がレーダー波をどれだけ反射するかを表す指標)に変換する必要があります。

# 1. Conversion to Power
caldB = -83
calPwr = np.power(10.,caldB/10.)
rasterstack_pwr = np.power(rasterstack,2.)*calPwr

# 2. Compute Means
rs_means_pwr = np.mean(rasterstack_pwr,axis=(1,2))

# 3. Convert to dB
rs_means_dB = 10.*np.log10(rs_means_pwr)

変換、計算の流れは以下です。

• パワー値への変換
  • caldB = -83: 校正値（dB単位）を設定します。SARデータの信号強度は、この値に基づいて補正されます。
  • calPwr = np.power(10., caldB / 10.): 校正値をパワー値に変換します。dB単位の値をパワー値に変換するには、10を底とする指数関数を使用します。
  • rasterstack_pwr = np.power(rasterstack, 2.) * calPwr: SARデータの各ピクセルの強度をパワー値に変換し、校正値で補正します。各ピクセルの値を2乗することで、信号の強度をパワー（電力）に変換します。
• 平均値の計算
  • np.mean(rasterstack_pwr, axis=(1,2)): 全ての画像のパワー値に対して、各バンド（チャンネル）の平均パワーを計算します。axis=(1,2) を指定することで、各バンドの全ピクセルの平均を求めます。
• dBへの変換
  • 10 * np.log10(rs_means_pwr): 平均パワー値をdB単位に変換します。平均パワー値を対数に変換し、その後10倍します。これにより、平均パワーがdB単位で表されます。

取得された平均値の数を表示します。
band毎に平均値を計算したため、tindexに格納されていた日付の数と同じ値が表示されます。

rs_means_pwr.shape # Check that we got the means over time

取得された平均値をグラフで表示します。

# 4. Plot
fig=plt.figure(figsize=(16,4))
ax1=fig.add_subplot(111)
ax1.plot(tindex,rs_means_pwr)
ax1.set_xlabel('Date')
ax1.set_ylabel('$\overline{\gamma^o}$ [power]')
ax2=ax1.twinx()
ax2.plot(tindex,rs_means_dB,color='red')
ax2.set_ylabel('$\overline{\gamma^o}$ [dB]')
fig.legend(['power','dB'],loc=1)
plt.title('Time series profile of average band backscatter $\gamma^o$ ')

日付ごとの平均値をリストで表示します。

a =pd.Series(rs_means_dB,index=tindex)

#This will print a list of the dates and the respective dB mean values
a

A two part figure with moving global mean backscatter of the time series in dB

上で表示したアニメーション画像と連動してグラフが更新されていくように表示します。
コードは上述のものを組み合わせたものなので、ここでは解説を割愛します。

%%capture
import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np

fig, (ax1,ax2) = plt.subplots(1,2,figsize=(16,4),gridspec_kw = {'width_ratios':[1, 3]})
vmin=np.percentile(rasterstack.flatten(),5)
vmax=np.percentile(rasterstack.flatten(),95)
im = ax1.imshow(raster0,cmap='gray',vmin=vmin,vmax=vmax)
ax1.set_title("{}".format(tindex[0].date()))
ax1.set_axis_off()
ax2.axis([tindex[0],tindex[-1],rs_means_dB.min(),rs_means_dB.max()])
ax2.set_ylabel('$\overline{\gamma^o}$ [dB]')
ax2.set_xlabel('Date')
ax2.set_ylim((-15,-5))
l, = ax2.plot([],[])

def animate(i):
    ax1.set_title("{}".format(tindex[i].date()))
    im.set_data(rasterstack[i])
    ax2.set_title("{}".format(tindex[i].date()))
    l.set_data(tindex[:(i+1)],rs_means_dB[:(i+1)])

# Interval is given in milliseconds
ani = matplotlib.animation.FuncAnimation(fig, animate, frames=rasterstack.shape[0], interval=400)

from IPython.display import HTML
HTML(ani.to_jshtml())

まとめ

本記事では「SAR Handbook」の3章のトレーニングPart2について写経し、自分なりに理解してみました。
どちらかというと自分のための備忘録として記載しましたが、SARについて勉強し始めたところですので、もし間違っている点等ありましたらコメントいただけると幸いです。