はじめに

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

以下の記事の続きで、本記事ではChapter3のPart4について記載します。
https://qiita.com/drafts/71ef7555f9497c9041f7/

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 4 - SAR Time Series Change Point Detection

In this chapter we introduce the advanced concepts of change point detection in time series. One of the goals of change detection for forest applications is to identifydisturbance over an observation period and the timing of events. Many tools for change point detection stem from the financial sector and are available today with differentcomplexities. In this workbook we will analyze time series signatures from SAR with emphasis on forest time series. We will start by exploring time series at pixel levels andwill work up to a change point detection scenario with image based analysis.

本章では、時系列における変化点検出を紹介します。
森林への応用における変化点検出の目的の1つは、観測期間中の撹乱とイベントのタイミングを特定することです。変化点検出のための多くのツールは金融セクターから派生したものであり、現在は様々なものが利用可能です。
このワークブックでは、森林時系列に重点を置いて、SARデータの時系列分析します。まず、ピクセルレベルの時系列を探索することから始め、画像ベースの解析による変化点検出のシナリオに取り組みます。

まずはライブラリをimportします。

# Importing relevant python packages
import pandas as pd
import gdal
import numpy as np
import time,os

# For plotting
%matplotlib inline
import matplotlib.pylab as plt
import matplotlib.patches as patches

font = {'family' : 'monospace',
          'weight' : 'bold',
          'size'   : 18}
plt.rc('font',**font)

Select the Project data set and time series data

データをダウンロードします。

West Africa - Niamey Deforestation Site

datadirectory='/Users/rmuench/Downloads/wa/cra/'
datefile='S32631X402380Y1491460sS1_A_vv_0001_A_mtfil.dates'
imagefile='S32631X402380Y1491460sS1_A_vv_0001_A_mtfil.vrt'

os.chdir(datapath)

Acquisition Dates

Read from the Dates file the dates in the time series and make a pandas date index

datesファイルから日付を読み込んで、pandasで日付のindexを作成します。

dates=open(datefile).readlines()
tindex=pd.DatetimeIndex(dates)
j=1
print('Bands and dates for',imagefile)
for i in tindex:
    print("{:4d} {}".format(j, i.date()),end=' ')
    j+=1
    if j%5==1: print()

Image data

Get the time series raster stack from the entire training data set.

スタックされたラスターファイル(複数の日付の画像)を配列で読み込みます。

rasterstack=gdal.Open(imagefile).ReadAsArray()

Data Pre-Processing

Plot the global means of the Time Series

Conversion to power

Compute means

Convert to dB

Make a pandas time series

Plot time series of means

平均値の時系列変化をプロットします。
手順は以下です。

パワー値への変換
平均値の計算
dB値への変換
pandasで時系列データの作成
平均値の時系列変化のプロット

まずは、1.~3.を行います。
part1でも記載がありましたが、平均値の計算はパワー値で行うという点に注意してください。dB値での平均値の計算は異なる結果となります。

# 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)

4.を実施します。

# 4. Make a pandas time series object
ts = pd.Series(rs_means_dB,index=tindex)

次に5.です。上で作成していた日付のindexと平均値を俺線グレフにプロットしていきます。

# 5. Use the pandas plot function of the time series object to plot
# Put band numbers as data point labels
plt.figure(figsize=(16,8))
ts.plot()
xl = plt.xlabel('Date')
yl = plt.ylabel('$\overline{\gamma^o}$ [dB]')
for xyb in zip(ts.index,rs_means_dB,range(1,len(ts)+1)):
    plt.annotate(xyb[2],xy=xyb[0:2])

zip(ts.index, rs_means_dB, range(1, len(ts) + 1))で、tsのインデックス、rs_means_dB、およびバンド番号のリストを組み合わせたタプルを生成します。
ts.indexは時系列データの日時インデックスです。
rs_means_dBは後方散乱係数の平均値（dB単位）のリストです。
range(1, len(ts) + 1)は1からデータポイント数までのバンド番号を生成します。
plt.annotate(xyb[2], xy=xyb[0:2])で、各データポイントにバンド番号をラベルとして追加します。
- xyb[2]はバンド番号です。
- xyb[0:2]はラベルの位置を示す座標（日時インデックスと平均値のペア）です。

EXERCISE

Look at the global means plot and determine from the tindex array at which dates you see maximum and minimum values. Are relative peaks associated with seasons?

「平均値のプロットを見て、tindex配列から、どの日付で最大値と最小値を見るかを決定します。相対的なピークは季節と関連しているか。」と問われています。
こちらは、回答/解説がないのですが、グラフを見ると2016年も2017年も7月から11月の間にピークがあり、季節と関連していると思われます。

Generate Time Series for Point Locations or Subsets

In python we can use the matrix slicing rules (Like Matlab) to obtain subsets of the data. For example to pick one pixel at a line/pixel location and obtain all band values, use:
[:,line,pixel] notation.
Or, if we are interested in a subset at an offset location we can use:
[:,yoffset:(yoffset+yrange),xoffset:(xoffset+xrange)]
In the section below we will learn how to generate time series plots for point locations (pixels) or areas (e.g. a 5x5 window region). To show individual bands, we define a showImage function which incorporates the matrix slicing from above.

行列のスライスを使ってデータのサブセットを取得します。例えば、あるライン/ピクセルの位置で1つのピクセルを選び、すべてのバンド値を得るには、次のようにします：
[:,line,pixel]
また、オフセット位置のサブセットを取得する場合は、次のようにします。
[:,yoffset:(yoffset+yrange),xoffset:(xoffset+xrange)]
以下のセクションでは、点の位置（ピクセル）またはエリア（例えば5x5のウィンドウ領域）の時系列プロットを生成する方法を学びます。個々のバンドを表示するために、上述のマトリックススライスを組み込んだshowImage関数を定義します。

def showImage(rasterstack,tindex,bandnbr,subset=None,vmin=None,vmax=None):
    '''Input: 
    rasterstack stack of images in SAR power units
    tindex time series date index
    bandnbr bandnumber of the rasterstack to dissplay'''
    fig = plt.figure(figsize=(16,8))
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)
    
    # If vmin or vmax are None we use percentiles as limits:
    if vmin==None: vmin=np.percentile(rasterstack[bandnbr-1].flatten(),5)
    if vmax==None: vmax=np.percentile(rasterstack[bandnbr-1].flatten(),95)

    ax1.imshow(rasterstack[bandnbr-1],cmap='gray',vmin=vmin,vmax=vmax)
    ax1.set_title('Image Band {} {}'.format(bandnbr,tindex[bandnbr-1].date()))
    if subset== None:
        bands,ydim,xdim=rasterstack.shape
        subset=(0,0,xdim,ydim)
        
    ax1.add_patch(patches.Rectangle((subset[0],subset[1]),subset[2],subset[3],fill=False,edgecolor='red'))
    ax1.xaxis.set_label_text('Pixel')
    ax1.yaxis.set_label_text('Line')
    
    ts_pwr=np.mean(rasterstack[:,subset[1]:(subset[1]+subset[3]),
                       subset[0]:(subset[0]+subset[2])],axis=(1,2))
    ts_dB=10.*np.log10(ts_pwr)
    ax2.plot(tindex,ts_dB)
    ax2.yaxis.set_label_text('$\gamma^o$ [dB]')
    ax2.set_title('$\gamma^o$ Backscatter Time Series')
    # Add a vertical line for the date where the image is displayed
    ax2.axvline(tindex[bandnbr-1],color='red')

    fig.autofmt_xdate()

この関数showImage()は、スタックされたSARデータから特定のバンドの画像を表示し、そのバンドの領域に対応する後方散乱値の時系列変化をプロットします。

関数の引数は以下です。
- rasterstack：SARパワーユニットで表された画像のスタック（3次元配列）。
- tindex：時系列の日付インデックス（pandas.DatetimeIndexのようなもの）。
- bandnbr：表示する画像スタックのバンド番号（1から始まる）。
- subset：画像のサブセットを指定するためのタプル（(x, y, width, height)）。デフォルトはNoneで、全体を表示。
- vminとvmax：表示する画像のピクセル値の最小値と最大値。デフォルトはNoneで、パーセンタイルで設定。
vminもしくはvmaxが指定されていない場合は、画像の5パーセントタイルと95パーセントタイルをそれぞれ最大値最小値に設定しています。
subsetが指定されていない場合は、画像全体をサブセットとして設定しています。subsetの領域は赤枠で囲うようになっています。
時系列データの計算配下の手順で行います。
- subset領域のパワー値の平均値を計算します。
- 算出した平均値をdB単位に変換します。

実際に使用してみます。

bandnbr=24
subset=[5,20,3,3]
showImage(rasterstack_pwr,tindex,bandnbr,subset)

bandnbr=43
showImage(rasterstack_pwr,tindex,bandnbr,subset)

bandnbr=43はbandnbr=24に比べて赤枠の領域が全体的に暗く見え、実際に平均値も低くなっています。

EXERCISE

For subset (5,20,3,3):
1.What are the dates where backscatter falls below -11 dB?
2.Compute the gradients (backscatter difference between two consecutive dates.
3.What is the largest gradient of backscatter drop between two consecutive dates?
4.What are the dates associated with this gradient (before and after)?

質問/課題は以下です。
・後方散乱が-11dBを下回った日は？
・勾配（連続する2つの日付間の後方散乱の差）を計算する。
・連続する2つの日付間の後方散乱の低下の最大の勾配は？
・この勾配に関連する日付（前後）は？

筆者はデータをダウンロードできていないため、実際に課題に取り組むことはできないのですが、以下のコードで１つ隣との差分を求めてグラフをプロットすることができます。こちらで後方散乱の低下が最大のところを見つけることができるでしょう。

gradient_lag1 = ts.diff(1)
gradient_lag1.plot()

Helper function the generate a time series object

def timeSeries(rasterstack_pwr,tindex,subset,ndv=0.):
# Extract the means along the time series axes
# raster shape is time steps, lines, pixels.
# With axis=1,2, we average lines and pixels for each time
# step (axis 0)
raster=rasterstack_pwr.copy()
if ndv != np.nan: raster[np.equal(raster,ndv)]=np.nan
ts_pwr=np.nanmean(raster[:,subset[1]:(subset[1]+subset[3]),
subset[0]:(subset[0]+subset[2])],axis=(1,2))
# convert the means to dB
ts_dB=10.*np.log10(ts_pwr)
# make the pandas time series object
ts = pd.Series(ts_dB,index=tindex)
# return it
return ts

この関数 timeSeries は、与えられたSARデータのスタックから特定のサブセット領域の時系列データを抽出し、そのデータをdB単位に変換してPandasの時系列オブジェクトとして返します。処理の流れは以下です。

無効値の処理
- rasterstack_pwrのコピーを作成してrasterに保存しています。
- ndvが指定されている場合、その値をnp.nanに置き換えます。これにより、無効値が計算に含まれないようにします。
サブセット領域の平均値を計算
- subset[1]:(subset[1]+subset[3])は行範囲。
- subset[0]:(subset[0]+subset[2])は列範囲。
- axis=(1,2)は、行と列の平均を取って、時間軸ごとに平均を計算します。
算出した平均値をdB単位に変換
Pandasのデータフレームを作成

実際に関数を使って時系列データを作成してみましょう。

ts = timeSeries(rasterstack_pwr,tindex,subset)
_=ts.plot(figsize=(16,4)) # _= is a trick to suppress more output.

Seasonal Subsets of time series records

Let's expand upon SAR time series analysis. Often it is desirable to subset time series by season or months to compare with similar conditions of a previous year's observation. For example, in analyzing C-Band backscatter data, it might be useful to limit comparative analysis to dry season observations only as soil moisture might confuse signals during the wet seasons. In this section we will expand upon the concepts of subsetting time series along the time axis. We will make use of the pandas datatime index tools:

Month
Day of year

First we extract the time series again for a area at the subset location (5,20,5,5). We then convert the pandas time series to a pandas DataFrame to allow for more processing options. We also label the data value column as 'g0' for gamma0:

SARの時系列解析について説明します。
多くの場合、季節や月ごとにサブセットとして、前年の同じような観測条件のデータと比較することが望ましいです。例えば、Cバンドの後方散乱データを解析する場合、雨季の間は土壌水分が信号を乱す可能性があるため、乾季の観測データのみに限定して比較解析することが有効な場合があります。
このセクションでは、時系列データを時間に沿ってサブセット化することを紹介します。まず、サブセット位置(5,20,5,5)のエリアの時系列データを再度抽出します。そして、pandasの時系列データをpandas DataFrameに変換し、より多くの処理オプションを可能にします。また、データ値の列には、gamma0 を表す 'g0' というラベルを付けます。

subset=(5,20,5,5)
ts = timeSeries(rasterstack_pwr,tindex,subset)
tsdf = pd.DataFrame(ts,index=ts.index,columns=['g0'])
# Plot
ylim=(-20,-5)
tsdf.plot(figsize=(16,4))
plt.title('Sentinel-1 C-VV Time Series Backscatter Profile, Subset: 5,20,5,5 ')
plt.ylabel('$\gamma^o$ [dB]')
plt.ylim(ylim)
_=plt.legend(["C-VV Time Series"])

Start the time series in November 2015

We can use the pandas index parameters like month to make seasonal subsets

monthのようなpandasのインデックスパラメータを使って、季節ごとのサブセットを作ることができます。

以下のコードでは2015年11月1日以降のデータを抽出してグラフを表示しています。

tsdf_sub1=tsdf[tsdf.index>'2015-11-01']

# Plot
tsdf_sub1.plot(figsize=(16,4))
plt.title('Sentinel-1 C-VV Time Series Backscatter Profile, Subset: {}'.format(subset))
plt.ylabel('$\gamma^o$ [dB]')
plt.ylim(ylim)
_=plt.legend(["C-VV Time Series"])

Subset by months

We can make use of pandas DateTimeIndex object index.month and numpy's logical_and function to subset a time series easily by month.

pandasのDateTimeIndexオブジェクトindex.monthとnumpyのlogical_and関数を利用することで、時系列を月ごとに簡単にサブセットすることができます。

まずは3月から5月のデータを抽出し、そのデータをプロットします。

tsdf_sub2=tsdf_sub1[
    np.logical_and(tsdf_sub1.index.month>=3,tsdf_sub1.index.month<=5)]

# Plot
fig, ax = plt.subplots(figsize=(16,4))
tsdf_sub2.plot(ax=ax)
plt.title('Sentinel-1 C-VV Time Series Backscatter Profile, Subset: {}'
          .format(subset))
plt.ylabel('$\gamma^o$ [dB]')
plt.ylim(ylim)
_=plt.legend(["March-May"])

次にその他の月のデータを抽出してみます。
numpyのブーリアン反転関数を使えば、選択範囲を反転させることができます。

tsdf_sub3=tsdf_sub1[np.invert(
    np.logical_and(tsdf_sub1.index.month>=3,tsdf_sub1.index.month<=5))]

# Plot
fig, ax = plt.subplots(figsize=(16,4))
tsdf_sub3.plot(ax=ax)
plt.title('Sentinel-1 C-VV Time Series Backscatter Profile, Subset: {}'
          .format(subset))
plt.ylabel('$\gamma^o$ [dB]')
plt.ylim(ylim)
_=plt.legend(["June-February"])

Group time series by year to compare average backscatter values

後方散乱値の平均値を比較するため、時系列データを年ごとにグループ化します。

ts_sub_by_year = tsdf_sub1.groupby(pd.Grouper(freq="Y"))

コードの詳細は以下です。

pd.Grouper(freq="Y")：PandasのGrouperオブジェクトを使用して、データを年ごとにグループ化します。
groupby：Grouperを使用して、DataFrameを年ごとにグループ化します。

年ごとにグループ化したデータをプロットします。

fig, ax = plt.subplots(figsize=(16,4))
for label, df in ts_sub_by_year:
df.g0.plot(ax=ax, label=label.year)
plt.legend()
# ts_sub_by_year.plot(ax=ax)
plt.title('Sentinel-1 C-VV Time Series Backscatter Profile, Subset: {}'
.format(subset))
plt.ylabel('$\gamma^o$ [dB]')
plt.ylim(ylim)

Make a pivot table to group year and sort by day of year for plotting overlapping time series

重なり合う時系列データをプロットするために、年をグループ化し、日ごとにソートするピボット・テーブルを作成する。

まず、データフレームに日と年の2つの列を追加します。

# Add doy
tsdf_sub1 = tsdf_sub1.assign(doy=tsdf_sub1.index.dayofyear)
# Add year
tsdf_sub1 = tsdf_sub1.assign(year=tsdf_sub1.index.year)

元の時系列データフレーム (tsdf_sub1) に「day of year (doy)」と「year」という2つの新しい列を追加しています。これにより、各データポイントがどの年の何日目に属するかがわかるようになります。

Add doy
- tsdf_sub1.index.dayofyear：データフレームのインデックス（日時）の各要素について、その年の1月1日からの経過日数を取得します。
- assign(doy=...)：取得した経過日数を新しい列「doy」としてデータフレームに追加します。
Add year
- tsdf_sub1.index.year：データフレームのインデックス（日時）の各要素について、その年を取得します。
- assign(year=...)：取得した年を新しい列「year」としてデータフレームに追加します。

さらに、曜日をインデックス、年をカラムとするピボット・テーブルを作成します。

piv=pd.pivot_table(tsdf_sub1,index=['doy'],columns=['year'],values=['g0'])
# Set the names for the column indices
piv.columns.set_names(['g0','Year'],inplace=True)
print(piv.head(10))
print('...\n',piv.tail(10))

カラム名を設定します。

piv.columns.set_names(['g0','year'],inplace=True)

As we can see, there are NaN (Not a Number) values on the days in a year where no acquisition took place. Now we use time weighted interpolation to fill the dates for all the observations in any given year. For time weighted interpolation to work we need to create a dummy year as a date index, perform the interpolation, and reset the index to the day of year. This is accomplished with the following steps:

取得が行われなかった年の日にNaN(Not a Number)値があリます。ここで、任意の年の全ての日付を埋めるために、時間加重補間を使用します。時間重み付き補間のために、日付インデックスとしてダミーの年を作成し、補間を実行し、インデックスをその年のその日にリセットする必要があります。

# Add fake dates for year 100 to enable time sensitive interpolation 
# of missing values in the pivot table
year_doy = ['2100-{}'.format(x) for x in piv.index]
y100_doy=pd.DatetimeIndex(pd.to_datetime(year_doy,format='%Y-%j'))

# make a copy of the piv table and add two columns
piv2=piv.copy()
piv2=piv2.assign(d100=y100_doy) # add the fake year dates
piv2=piv2.assign(doy=piv2.index) # add doy as a column to replace as index later again

# Set the index to the dummy year
piv2.set_index('d100',inplace=True,drop=True)

# PERFORM THE TIME WEIGHTED INTERPOLATION
piv2 = piv2.interpolate(method='time')  # TIME WEIGHTED INTERPOLATION!

# Set the index back to day of year.
piv2.set_index('doy',inplace=True,drop=True)

詳細な流れは以下です。

架空の日付を追加
- year_doy：pivテーブルのインデックス（日数）を使って、2100年の日付リストを作成します。
- pd.to_datetime(year_doy, format='%Y-%j')：リスト内の日付文字列をDatetimeIndexに変換します。'%Y-%j'は、年とその年の日数を表すフォーマットです。
- y100_doy：これにより、2100年の日付のインデックスが生成されます。
ピボットテーブルのコピーを作成し、列を追加
- piv.copy()：pivテーブルのコピーを作成します。
- piv2.assign(d100=y100_doy)：d100列に2100年の日付を追加します。
- piv2.assign(doy=piv2.index)：doy列に元のインデックス（day of year）を追加します。
インデックスを架空の年に設定
- piv2.set_index('d100', inplace=True, drop=True)：d100列を新しいインデックスとして設定します。この際、元のd100列は削除されます。
時間を考慮した補間の実行
- piv2.interpolate(method='time')：時間を考慮した補間を実行します。この方法では、時間軸に沿ったデータポイントを補間します。
インデックスを元のday of yearに戻す
- piv2.set_index('doy', inplace=True, drop=True)：元のdoy（day of year）をインデックスに戻します。この際、doy列は削除されます。

Let's inspect the new pivot table and see wheather we interpolated the NaN values where it made sense:

新しいピボットテーブルを確認し、NaN値の補間がされたかみてみましょう。

print(piv2.head(10))
print('...\n',piv2.tail(10))

piv2.plot(figsize=(16,8))
plt.title('Sentinel-1 C-VV Time Series Backscatter Profile,\
Subset: 5,20,5,5 ')
plt.ylabel('$\gamma^o$ [dB]')
plt.xlabel('Day of Year')
_=plt.ylim(ylim)

Change Detection on the Time Series Data

We can now analyze the time series for change. We will discuss two approaches:

Year-to-year differencing of the subsetted time series

Cumulative Sum based change detection

時系列の変化を分析します。ここでは2つのアプローチについて説明します。
1.サブセット化した時系列の年別差分分析
2.累積和に基づく変化検出

Year-to-Year Change Detection

We compute the differences between the interpolated time series and look for change with a threshold value.
補間された時系列間の差分を計算し、閾値による変化を検出します。

差分の分析に使用する閾値を設定します。
今回閾値には3dBを設定していますが、3dBの変化はSARの反射が2倍または半分になっていることを意味します。

# Difference between years
# Set a dB change threshold
thres=3

差分の計算を実施し、グラフを表示します。

diff1716 = (piv2.g0[2017]-piv2.g0[2016])
_=diff1716.plot('line')

変化の絶対値が閾値以上のdoyを抽出します。

thres_exceeded = diff1716[abs(diff1716) > thres]
thres_exceeded

From the three_exceeded dataframe we can infer the first date at which the threshold was exeeded. We would label that as a change point. As an additional criteria for labeling a change point, one can also consider the number of observations after identification of a change point where backscatter differed from the year before. If only one or two observations differed from the year before this could be considered an outlier. Addtionally, one can introduce smoothing operations with the interpolation

three_exceededから、閾値を超えた最初の日付を推測し、それを変化点としてラベル付けします。変化点をラベル付けするための追加基準として、後方散乱が前年と異なる変化点の数を考慮することも考えられます。もし1つか2つの日付だけが前年と異なるなら、これは異常値とみなされるかもしれません。

Cumulative Sums for Change Detection

Another approach to detect change in regularly acquired data is employing cumulative sums. Changes are determined against mean observations of time series. A full explanation and examples from the the financial sector can be found at http://www.variation.com/cpa/tech/changepoint.html

定期的に取得されるデータの変化を検出するもう一つのアプローチは、累積和を採用することです。変化は時系列の平均観測値に対して決定されます。詳細な説明と金融業界での例は、http://www.variation.com/cpa/tech/changepoint.html で紹介されています。

Time Series and means

First let's consider a time series and it's mean observation. We look at two full years of observations from Sentinel-1 data for an area where we suspect change. In the following we consider 𝑋 as a time series

まず、時系列データの平均値を考えてみます。ここでは、変化が疑われる地域のSentinel-1データから2年分のデータを見ます。以下では、𝑋を時系列として考えます。

$$X = (X_1, X_2, \ldots, X_n)$$
with

$X_i$ SAR backscatter at time $i = 1, \ldots, n$
$n$ number of observations in the time series

subset=(5,20,3,3)
#subset=(12,5,3,3)
ts1 = timeSeries(rasterstack_pwr,tindex,subset)
X = ts1[ts1.index>'2015-10-31']

Filtering the time series for outliers

It is advantageous in noisy SAR time series data like C-Band data to filter on the time axis. Pandas offers a "rolling" function for these purposes. With that function we can choose, for example, a median filter along the time axis. Below is an example of a median filter for an observation filters the time series when the observation before and after a time stamps are part of the filter.

Cバンドのようなノイズの多いSARデータでは、時間軸でフィルタリングするのが良いです。Pandasはこのような目的のためにrolling関数を提供しています。この機能により、例えば時間軸に沿った中央値フィルタを選択することができます。以下は、ある観測値に対する中央値フィルタの例で、あるタイムスタンプの前後の観測値がフィルタの一部となるように時系列をフィルタリングします。

Xr=X.rolling(5,center=True).median()
Xr.plot()
_=X.plot()

Let's plot the time series and it's mean over the time span

時系列平均をプロットしてみます。

X=Xr  # Uncomment if rolling mean is wanted for further computation
Xmean = X.mean()
fig,ax=plt.subplots(figsize=(16,4))
X.plot()
plt.ylabel('$\gamma^o$ [dB]')
ax.axhline(Xmean,color='red')
_=plt.legend(['$\gamma^o$','$\overline{\gamma^o}$'])

上で作成したローリングメディアンの時系列データのプロットを作成し、データの平均値を赤い水平線として表示しています。

Let's determine the residuals of the time series against the mean

平均に対する時系列データの残差を求めます。
$$R = X_i - \overline{X}$$

R = X - Xmean

Now we compute the cumulative sum of the residuals and plot it

ここで残差の累積和を計算し、プロットします。
$$S = \sum_{n=1}^{n}R_i$$

S = R.cumsum()
_=S.plot(figsize=(16,8))

An estimator for the magnitude of change is given as the difference between the maximum and minimum value of S

変化の大きさの推定値は、Sの最大値と最小値の差として求められます。
$$S_{DIFF} = S_{MAX} - S_{MIN}$$

Sdiff=S.max() - S.min()
Sdiff

A candidate change point is identified from the 𝑆 curve at the time where 𝑆𝑀𝐴𝑋 is found:

変化点の候補は、$S_{MAX}$が見つかった時点の𝑆曲線から特定される。
$$T_{CP_{before}} = T(S_i = S_{MAX})$$
with
*$T_{CP_{before}}$ Timestamp of last observation before change
*$S_i$ Cumulative Sum of R with $i = (1, \ldots, n)$

𝑛 Number of observations in the time series

変化後の最初の観測点$T_{CP_{before}}$は、$T_{CP_{before}}$に続く時系列の最初の観測点として求められる。
時系列Xの例では、以下のようになる。

t_cp_before = S[S==S.max()].index[0]
print('Last date before change: {}'.format(t_cp_before.date()))

t_cp_after = S[S.index > t_cp_before].index[0]
print('First date after change: {}'.format(t_cp_after.date()))

Bootstrapping the cumulative sums by randomly reordering the time series

We can determine if an identified change point is indeed a valid detection by randomly reordring the time series and comparing the various S curves. During bootstrapping we count how many times the $S_{DIFF}$ values are greater than $S_{DIFF_{rondom}}$ of the identified change point. A confindence level 𝐶𝐿 is computed as

時系列をランダムに再検索し、様々なSカーブを比較することで、特定された変化点が本当に有効な検出かどうかを判断することができます。ブートストラップ中に、$S_{DIFF_{rondom}}$値が識別された変化点の$S_{DIFF}$より大きい回数を数えます。確信度𝐶𝐿は次のように計算されます。

$CL = \frac{N_{GT}}{N_{bootstraps}}$
with

$N_{GT}$ Number of times $S_{DIFF} > S_{DIFF_{rondom}}$
$N_{bootstraps}$ Number of bootstraps randomizing R

Another metric for the significance of a change point is 1 minus the ratio of the mean of the 𝑆𝐷𝐼𝐹𝐹𝑟𝑎𝑛𝑑𝑜𝑚 values and 𝑆𝐷𝐼𝐹𝐹 . The closer this value is approaching 1, the more significant the change point:

重要な変化点を示すもう1つの指標は、1から$S_{DIFF_{rondom}}$の平均値の比率を引いた値です。この値が1に近ければ近いほど、その変化点は有意です。

The python code to conduct the boot strapping, including visualization of the S curves is below:

S字曲線の視覚化を含むブートストラッピングを行うためのコードは以下の通りです。

n_bootstraps=500  # bootstrap sample size
fig,ax = plt.subplots(figsize=(16,8))
S.plot(ax=ax,linewidth=3)
ax.set_ylabel('Cumulative Sums of the Residuals')
fig.legend(['S Curve for Candidate Change Point'],loc=3)
Sdiff_random_sum=0
Sdiff_random_max=0  # to keep track of the maxium Sdiff of the 
               # bootstrapped sample
n_Sdiff_gt_Sdiff_random=0  # to keep track of the maxium Sdiff of the 
               # bootstrapped sample
for i in range(n_bootstraps):
    Rrandom = R.sample(frac=1)  # Randomize the time steps of the residuals
    Srandom = Rrandom.cumsum()
    Sdiff_random=Srandom.max()-Srandom.min()
    Sdiff_random_sum += Sdiff_random
    if Sdiff_random > Sdiff_random_max:
        Sdiff_random_max = Sdiff_random
    if Sdiff > Sdiff_random:
        n_Sdiff_gt_Sdiff_random += 1
    Srandom.plot(ax=ax)
_=ax.axhline(Sdiff_random_sum/n_bootstraps)

ここのコードは少し難しいと感じたので、内容を詳細に解説したいと思います。

1.初期設定とプロットの準備

n_bootstraps = 500  # ブートストラップサンプルのサイズ
fig, ax = plt.subplots(figsize=(16,8))
S.plot(ax=ax, linewidth=3)
ax.set_ylabel('Cumulative Sums of the Residuals')
fig.legend(['S Curve for Candidate Change Point'], loc=3)

n_bootstrapsはブートストラップサンプルのサイズ（500）を指定しています。
fig, axはプロットを作成するためのFigureとAxesオブジェクトです。
S（累積和時間シリーズオブジェクト）をプロットし、y軸を設定します。
fig.legendで凡例を設定します。

2.ブートストラップ法の実装

Sdiff_random_sum = 0
Sdiff_random_max = 0  # ブートストラップサンプルの最大Sdiffを記録するための変数
n_Sdiff_gt_Sdiff_random = 0  # ブートストラップサンプルの最大Sdiffを記録するための変数

for i in range(n_bootstraps):
    Rrandom = R.sample(frac=1)  # 残差のタイムステップをランダム化
    Srandom = Rrandom.cumsum()
    Sdiff_random = Srandom.max() - Srandom.min()
    Sdiff_random_sum += Sdiff_random
    if Sdiff_random > Sdiff_random_max:
        Sdiff_random_max = Sdiff_random
    if Sdiff > Sdiff_random:
        n_Sdiff_gt_Sdiff_random += 1
    Srandom.plot(ax=ax)

Sdiff_random_sum、Sdiff_random_max、n_Sdiff_gt_Sdiff_randomは、それぞれ累積和の差の和、最大値、元の累積和の差を超える回数を記録するための変数です。
forループで 500 回のブートストラップサンプルを生成し、それぞれ累積してプロットします。
R.sample(frac=1)は残差データをランダムに並べ替えます。
Srandom = Rrandom.cumsum()でランダム化された残差の累積を意味します。
Sdiff_random = Srandom.max() - Srandom.min()でランダムな累積和の差が適用されます。
Sdiff_random_sumにこの差を追加し、最大値を変更します。
元の累積和の差 ( Sdiff) と比較し、大きくすればカウントします。

3.結果のプロット

_ = ax.axhline(Sdiff_random_sum/n_bootstraps)

ax.axhline で、ブートストラップサンプルの平均累積和の差を水平線としてプロットします。

検出された変化点の信頼レベル（Confidence Level, CL）を計算します。

n_Sdiff_gt_Sdiff_randomはブートストラップによって生成されたランダムな累積和の最大差（Sdiff_random）が、元の累積和の最大差（Sdiff）より大きかった回数をカウントしたものです。
n_Sdiff_gt_Sdiff_random / n_bootstraps: ブートストラップサンプルのうち、Sdiffより大きなSdiff_randomの割合を計算しています。この割合は、変化点がランダムである確率を示します。

CL = 1.*n_Sdiff_gt_Sdiff_random/n_bootstraps
print('Confidence Level for change point {} percent'.format(CL*100.))

今回はブートストラップによって生成されたランダムな累積和が、元の累和よりも全て大きかったため、信頼レベルは100%と算出されました。

次に、変化点の有意性指標（Change Point Significance Metric, CP_significance）を計算します。

有意性指標の計算は以下のように行います。
- 平均Sdiff_randomの計算: (Sdiff_random_sum / n_bootstraps)で、ブートストラップサンプルの平均Sdiff_randomを計算します。
- 元データのSdiffとブートストラップ平均Sdiff_randomの比率: (Sdiff_random_sum / n_bootstraps) / Sdiffは、元データの変化点の大きさと、ランダムサンプルの変化点の大きさの比率を示します。
- 有意性指標の計算: 1. - (Sdiff_random_sum / n_bootstraps) / Sdiffで、変化点の有意性指標を計算します。この指標は、1に近いほど変化点が有意であることを示します。

CP_significance = 1. - (Sdiff_random_sum/n_bootstraps)/Sdiff
print('Change point significance metric: {}'.format(CP_significance))

今回は0.591と算出されました。

Another useful metric to determine strength of a change point is the normalized integral $S_{ni}$　of the absolute values of the S curve:

変化点の強さを決定するもう一つの有用な指標は、S曲線の絶対値の正規化積分$S_{ni}$です。
$S_{normintegral} = \frac{\int_{i=1}^n \frac{abs(S_i)}{\max{abs(S)}} } {n}$

# NaN's to be excluded in the computation
S_ni=(S.abs()/S.abs().max()).cumsum().max()/len(S[S != np.nan])
print('Normalized Integral of cumulative sum: {}'.format(S_ni))

Selection of threshold values

$CL$ and $CP_{significance}$ can be used as threshold values for the acceptance or rejection of a candidate threshold. These values are to some degree specific to a SAR sensor and environmental conditions. E.g. L-Band SAR has a more pronounced decrease in backscatter after forest disturbance and logging, whereas C-Band can have more ambigious signals. Also moisture regime changes, e.g. with snow cover, freeze/thaw conditions or dry/wet season changes have an influence on the time series signal. For example El Nino years can suggest changes solely due to different wetting and dryup conditions pertinent to a particular year. For this reason other techniques can be added to the SAR time series ananlysis. Two techniques can readily be thought of:

Subsetting of time series by seasons

Detrending time series with global image means
If year-to-year comparison is the focus, the first approach likely leads to subsets that are too small for meaningful cumulative sum change point detection. The approach of interannual differencing as discussed above likely performs better.
In the following we explore the approach to detrend the data with global image means.

$CL$と$CP_{significance}$は、閾値候補の採否を決定するための閾値として用いることができる。これらの値はSARセンサーや環境条件にある程度依存します。例えば、LバンドSARは森林攪乱や伐採後の後方散乱の減少が顕著であるのに対し、Cバンドはより曖昧な信号を持つことがあります。また、積雪や凍結・融解、乾季・湿季の変化など、水分の変化も時系列信号に影響を与える。例えば、エルニーニョの年は、特定の年に関連する湿潤と乾燥の条件の違いだけによる変化を示唆することがある。このため、SAR時系列解析に他の手法を加えることができる。以下の2つの手法が考えられる：

季節による時系列のサブセット化
全球画像平均を用いた時系列のデトレンド処理

年ごとの比較に重点を置く場合、最初の方法ではサブセットが小さすぎて意味のある累積和の変化点を検出できない可能性が高い。上述したような年次間差分によるアプローチの方が、より良い結果が得られる可能性が高い。
以下では、大域的な画像平均を用いてデータをデトレンドするアプローチを検討する。

De-trending time series with global image means

The idea of de-trending time series with global image means should prepare time series for a somewhat more robust change point detection as global image time series anomalies stemming calibration or seasonal trends are removed prior to time series analysis. This de-trending needs to be performed with large subsets so real change is not influencing the image statistics.

NOTE: For our small subset, we will see some of these effects.

Let's start by building a global image means time series:

グローバル画像平均を用いて時系列をトレンド除去する手法は、いくらかロバストな変化点検出のために、校正や季節のトレンド起因するグローバル画像時系列の異常が時系列解析の前に除去として、準備する必要がある。このトレンド除去は、実際の変化が画像統計に影響を与えないように、大きなサブセットで実行する必要がある。
注：我々の小さなサブセットでは、これらの影響のいくつかが見られるだろう。

まず、グローバルな画像の時系列データを構築することから始める。

means_pwr = np.mean(rasterstack_pwr,axis=(1,2))
means_dB = 10.*np.log10(means_pwr)
gm_ts = pd.Series(means_dB,index=tindex)
gm_ts=gm_ts[gm_ts.index > '2015-10-31'] # filter dates
gm_ts=gm_ts.rolling(5,center=True).median()

このコードは時系列データの平均値をdBスケールに変換し、ローリングメディアンを計算するコードです。詳細は以下です。

パワー画像の平均値を計算
- rasterstack_pwrはSARデータのパワー画像の3次元配列です（時間、行、列の順）。
- 各時間ステップに対して、全ピクセルの平均パワー値を計算します。
- 結果は、各時間ステップに対応する1次元配列means_pwrになります。
平均値をdBスケールに変換
- パワー値をデシベル（dB）スケールに変換します。これは一般的な手法で、対数スケールを使うことで広範な値の範囲を扱いやすくします。
タイムシリーズデータとして格納
- dBスケールのデータをpandas.Seriesとして格納します。
- インデックスtindexは各時間ステップに対応する日付や時間の情報です。
特定の日付以降のデータにフィルタリング
- 2015年10月31日以降のデータにフィルタリングします。
ローリングメディアンを計算
- ウィンドウサイズ5でローリングメディアンを計算します。
- center=Trueにより、計算されたメディアン値がウィンドウの中央に配置されます。

gm_ts.plot()

X.plot()

Xd=X-gm_ts
Xmean=Xd.mean()
Xd.plot()

R = Xd - Xmean

Now we compute the cumulative sum of the residuals and plot it:

ここで残差の累積和を計算し、プロットします。
$$S = \sum_{n=1}^{n}R_i$$

S = R.cumsum()
_=S.plot(figsize=(16,8))

変化の大きさの推定値は、Sの最大値と最小値の差として求められます。
$$S_{DIFF} = S_{MAX} - S_{MIN}$$

Sdiff=S.max() - S.min()
Sdiff

変化点の候補は、$S_{MAX}$が見つかった時点の𝑆曲線から特定される。
$$T_{CP_{before}} = T(S_i = S_{MAX})$$
with
*$T_{CP_{before}}$ Timestamp of last observation before change
*$S_i$ Cumulative Sum of R with $i = (1, \ldots, n)$

𝑛 Number of observations in the time series

変化後の最初の観測点$T_{CP_{before}}$は、$T_{CP_{before}}$に続く時系列の最初の観測点として求められる。
時系列Xの例では、以下のようになる。

t_cp_before = S[S==S.max()].index[0]
print('Last date before change: {}'.format(t_cp_before.date()))

t_cp_after = S[S.index > t_cp_before].index[0]
print('First date after change: {}'.format(t_cp_after.date()))

Cumulative Sum Change Detection for the entire image

With numpy arrays we can apply the concept of cumulative sum change detection analysis effectively on the entire image stack. We take advantage of array slicing and axis-based computing in numpy. Axis 0 is the time domain in our raster stacks

numpyの配列を使えば、画像スタック全体に対して累積和変化検出分析の概念を効果的に適用することができます。numpyの配列スライスと軸ベースの計算を利用します。軸0はラスタスタックの時間領域です。

# Can do this in power or dB scale
X = rasterstack_pwr
# Filter out the first layer ( Dates >= '2015-11-1')
X_sub=X[1:,:,:]
tindex_sub=tindex[1:]
X= 10.*np.log10(X_sub)  # Uncomment to test dB scale

ここでは時系列データのうち、最初のデータ(band1)以外のデータを抽出し、dBスケールに変換しています。

plt.figure()
#Indicate the band number
bandnbr=0
vmin=np.percentile(X[bandnbr],5)
vmax=np.percentile(X[bandnbr],95)
plt.title('Band  {} {}'.format(bandnbr+1,tindex_sub[bandnbr].date()))
plt.imshow(X[0],cmap='gray',vmin=vmin,vmax=vmax)
_=plt.colorbar()

まずは最初のデータ(band1)を表示しました。

Xmean=np.mean(X,axis=0)
plt.figure()
plt.imshow(Xmean,cmap='gray')

次に上で作成した最初のデータ以外の平均値画像を表示しています。

R=X-Xmean

#Create an image that spatially displays the residuals (R)
plt.imshow(R[0])
plt.title('Residuals')
_=plt.colorbar()

これまで同様に時系列の各データから平均値画像を引いた残差を求めて、コンター図を作成しています。

S = np.cumsum(R,axis=0)
Smax= np.max(S,axis=0)
Smin= np.min(S,axis=0)
Sdiff=Smax-Smin
fig,ax=plt.subplots(1,3,figsize=(16,4))
vmin=Smin.min()
vmax=Smax.max()
p=ax[0].imshow(Smax,vmin=vmin,vmax=vmax)
ax[0].set_title('$S_{max}$')
ax[1].imshow(Smin,vmin=vmin,vmax=vmax)
ax[1].set_title('$S_{min}$')
ax[2].imshow(Sdiff,vmin=vmin,vmax=vmax)
ax[2].set_title('$S_{diff}$')
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
_=fig.colorbar(p,cax=cbar_ax)

累積和の最大値、最小値、およびその差を計算し、それらを画像として表示しています。

Mask Sdiff with a priori threshold for expected change

If we have an assumption as to how much actual change we expect in the image, we can threshold $S_{diff}$ to reduce computation of the bootstrapping. For land cover change we would not expect more than 5-10% change in a landscape. So, if the test region is reasonably large, setting a threshold for expected change to 10% would be appropriate. Thus we can set a mask with the 90th percentile of the histogram of $S_{diff}$. In our example we'll start out with a very conservative threshold of 50%.

The histogram for $S_{diff}$ is shown below:

もし、画像にどれくらいの実際の変化が予想されるかという仮定があれば、ブートストラップの計算を減らすために閾値Sを設定することができる。土地被覆の変化については、景観において5～10％以上の変化は期待できないと思われる。したがって、テスト領域がそれなりに大きければ、予想される変化の閾値を10%に設定するのが適切であろう。したがって、Sのヒストグラムの90パーセンタイルでマスクを設定する。今回の例では閾値を50パーセンタイルで設定する。
Sのヒストグラムを以下に示す。

#Display the Sdiff histogram
precentile=50
fig,ax=plt.subplots()
h=ax.hist(Sdiff.flatten(),bins=50)
thres=np.percentile(h[1],50)
print('At the {}% percentile, the threshold value is {:2.2f}'.format(precentile,
thres))
_=ax.axvline(thres,color='red')

累積和の差$S_{diff}$のヒストグラムを表示し、指定したパーセンタイル値に基づいて閾値を計算し、その閾値をヒストグラム上に赤い垂直線で表示しています。詳細を以下に示します。

ヒストグラムの表示
- Sdiff.flatten() は、Sdiff 配列を1次元に変換します。これにより、ヒストグラムを計算しやすくなります
- ax.hist(Sdiff.flatten(), bins=50) は、50個のビンを持つヒストグラムを作成します。h にはヒストグラムのデータが格納されます。
パーセンタイル値の計算
- np.percentile(h[1], 50) は、ヒストグラムのビン境界値 h[1] に対して50パーセンタイル（中央値）の値を計算します。
閾値をヒストグラムに表示
- ax.axvline(thres, color='red') は、計算された閾値を示す垂直線をヒストグラム上に赤色で描画します。

At the 50% percentile, the threshold value is ____ (printed above the histogram)

Using this threshold, we can visualize the candidate changepoints:

50％パーセンタイルでは、しきい値は19.82となりました。
この閾値を使って、変化点の候補を視覚化することができます。

Sdiffmask=Sdiff<thres
_=plt.imshow(Sdiffmask,cmap='gray')

求めた閾値を使用したマスク画像を表示しています。

Sdiffmask=Sdiff<thresで閾値以下をtrue、閾値以上をfalseとする配列を作成しています。
グレースケールで表示しているため、trueのピクセルが白、falseのピクセルが黒と表示されます。
よって、画像内の黒の部分が閾値以上のピクセルです。

Now we can filter our Residuals and perform bootstrapping analysis on these data. We make use of numpy masked arrays for this purpose.
残差をフィルタリングし、これらのデータに対してブートストラップ分析を実行します。numpyのマスク配列を使用します。

Rmask = np.broadcast_to(Sdiffmask,R.shape)
Rmasked = np.ma.array(R,mask=Rmask)

Sdiffmask が示す条件に基づいて R を部分的にマスクしています。

np.broadcast_to関数を使って、Sdiffmask を R の形状にブロードキャストします。これにより、Sdiffmask の形状が R と一致するように拡張されます。
np.ma.array を使用して、R をマスク配列 (np.ma.MaskedArray) に変換します。マスクは Rmask によって与えられます。つまり、Rmask が True の場所（つまり Sdiffmask が True の場所）はマスクされ、Rmasked の該当する要素は非表示にされます。

On the masked time series stack of residuals we can compute the cumulative sums

マスクされた残差の時系列スタックで累積和を計算します。

Smasked = np.ma.cumsum(Rmasked,axis=0)

np.ma.cumsum() 関数は、NumPy の masked array (np.ma.MaskedArray) に対して、指定した軸方向で累積和を計算するための関数です。この関数を使用すると、マスクされた要素を無視して累積和を計算することができます。

$S_{MAX}$, $S_{MIN}$, $S_{DIFF}$can also be computed on the masked arrays :

$S_{MAX}$, $S_{MIN}$, $S_{DIFF}$もマスクされた配列に対して計算します。

Smasked = np.ma.cumsum(Rmasked,axis=0)
Smasked_max= np.ma.max(Smasked,axis=0)
Smasked_min= np.ma.min(Smasked,axis=0)
Smasked_diff=Smasked_max-Smasked_min
fig,ax=plt.subplots(1,3,figsize=(16,4))
vmin=Smasked_min.min()
vmax=Smasked_max.max()
p=ax[0].imshow(Smasked_max,vmin=vmin,vmax=vmax)
ax[0].set_title('$S_{max}$')
ax[1].imshow(Smasked_min,vmin=vmin,vmax=vmax)
ax[1].set_title('$S_{min}$')
ax[2].imshow(Smasked_diff,vmin=vmin,vmax=vmax)
ax[2].set_title('$S_{diff}$')
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
_=fig.colorbar(p,cax=cbar_ax)

Bootstrapping over the masked change point candidates

We can now perform the bootstrapping analysis over the not masked out values. For efficient computing we permutate the index of the time axis

マスクされていない値に対してブートストラップ分析を行う。効率的な計算のために、時間s軸のインデックスを並べ替える。

random_index=np.random.permutation(Rmasked.shape[0])
Rrandom=Rmasked[random_index,:,:]
fig,ax=plt.subplots(1,2,figsize=(8,4))
ax[0].imshow(Rmasked[0])
ax[0].set_title('Band 0')
ax[1].imshow(Rrandom[0])
_=ax[1].set_title('Band 0 Randomized')

Rmaskedの時系列の順序をランダムにシャッフルした Rrandomを作成し、オリジナルの Rmasked の最初のバンドとランダム化されたRrandom の最初のバンドを並べて表示しています。

ランダムインデックスの生成
- np.random.permutation() は、指定したサイズの配列をランダムにシャッフルしたインデックスを生成します。ここでは、Rmasked の時間軸のサイズに基づいてランダムな順序を生成しています。
ランダム化された配列の作成
- ランダムなインデックスを使用して Rmasked の時間軸の順序をシャッフルし、新しい配列 Rrandom を作成しています。
プロットの表示

Smasked_max=np.ma.max(Smasked,axis=0)

Below is the numpy based implementation of the bootstrapping over all pixels. Note the efficient implementation using nympy masked arrays.

以下は、全ピクセルに対するブートストラップのnumpyベースの実装です。numpyのマスクされた配列を使った効率的な実装に注意してください。

n_bootstraps=1000  # bootstrap sample size

# to keep track of the maxium Sdiff of the bootstrapped sample:
Sdiff_random_max = np.ma.copy(Smasked_diff)
Sdiff_random_max[~Sdiff_random_max.mask]=0
# to compute the Sdiff sums of the bootstrapped sample:
Sdiff_random_sum = np.ma.copy(Smasked_diff)
Sdiff_random_sum[~Sdiff_random_max.mask]=0
# to keep track of the count of the bootstrapped sample
n_Sdiff_gt_Sdiff_random = np.ma.copy(Smasked_diff)
n_Sdiff_gt_Sdiff_random[~n_Sdiff_gt_Sdiff_random.mask]=0
for i in range(n_bootstraps):
    # For efficiency, we shuffle the time axis index and use that
    #to randomize the masked array
    random_index=np.random.permutation(Rmasked.shape[0])
    # Randomize the time step of the residuals
    Rrandom = Rmasked[random_index,:,:]
    Srandom = np.ma.cumsum(Rrandom,axis=0)
    Srandom_max=np.ma.max(Srandom,axis=0)
    Srandom_min=np.ma.min(Srandom,axis=0)
    Sdiff_random=Srandom_max-Srandom_min
    Sdiff_random_sum += Sdiff_random
    Sdiff_random_max[np.ma.greater(Sdiff_random,Sdiff_random_max)]=\
    Sdiff_random[np.ma.greater(Sdiff_random,Sdiff_random_max)]
    n_Sdiff_gt_Sdiff_random[np.ma.greater(Smasked_diff,Sdiff_random)] += 1

ブートストラップ法を用いてRmaskedの時間ステップをシャッフルし、累積和や最大値、最小値を計算し、それらの差（Sdiff_random）を元のSmasked_diffと比較して最大値や合計値のカウントを行っています。

初期化
- Sdiff_random_max、Sdiff_random_sum、および n_Sdiff_gt_Sdiff_random の初期化を行います。これらはそれぞれブートストラップサンプルの最大Sdiff、合計Sdiff、および元のSdiffがブートストラップサンプルのSdiffより大きい場合のカウントを追跡します。
ブートストラップサンプルの作成と計算
- 各ブートストラップサンプルについて、時間軸のインデックスをシャッフルして Rmasked をランダム化します。
- ランダム化された配列の累積和を計算し、その最大値と最小値を取得して Sdiff_random を計算します。
- Sdiff_random_sum に Sdiff_random を加算し、Sdiff_random_max を更新し、Smasked_diff が Sdiff_random より大きい場合のカウントを更新します。

Now we can compute for all pixels the confidence level CL, the change point significance metric $CP_{significance}$ and the product of the two as our confidence metric for identified changepoints

これですべてのピクセルの信頼度$CL$、変化点の有意性メトリック$CP_{significance}$、そして識別された変化点の信頼メトリックとしての2つの積を計算する。

CL = n_Sdiff_gt_Sdiff_random/n_bootstraps
CP_significance = 1.- (Sdiff_random_sum/n_bootstraps)/Sdiff
#Plot
fig,ax=plt.subplots(1,3,figsize=(16,4))
a = ax[0].imshow(CL*100)
fig.colorbar(a,ax=ax[0])
ax[0].set_title('Confidence Level %')
a = ax[1].imshow(CP_significance)
fig.colorbar(a,ax=ax[1])
ax[1].set_title('Significance')
a = ax[2].imshow(CL*CP_significance)
fig.colorbar(a,ax=ax[2])
_=ax[2].set_title('CL x S')

Now if we were to set a threshold of 0.5 for the product as identified change our change map would look like the following figure

ここで、もし閾値を0.5とし、変化として識別した場合、変更マップは以下の図のようになる。

cp_thres=0.5
plt.imshow(CL*CP_significance < cp_thres,cmap='cool')

Our last step is the idenficiaton of the change points is to extract the timing of the change. We will produce a raster layer that shows the band number of this first date after detected change. We will make use of the numpy indexing scheme. First, we create a combined mask of the first threshold and the identified change points after the bootstrapping. For this we use the numpy "mask_or" operation.

最後に、変化のタイミングを抽出する。検出された変化後の最初の日付のバンド番号を示すラスターレイヤーを作成する。numpyのインデックス・スキームを利用する。まず、最初の閾値と、ブートストラップ後に識別された変化点を組み合わせたマスクを作成する。これにはnumpyの "mask_or "演算を使用する。

# make a mask of our change points from the new threshold and the previous mask
cp_mask=np.ma.mask_or(CL*CP_significance<cp_thres,CL.mask)
# Broadcast the mask to the shape of the masked S curves
cp_mask2 = np.broadcast_to(cp_mask,Smasked.shape)
# Make a numpy masked array with this mask
CPraster = np.ma.array(Smasked.data,mask=cp_mask2)

新しい閾値（cp_thres）に基づいて、変更点のマスクを作成し、それをブロードキャストして累積和のマスク配列Smaskedに適用しています。最終的にこのマスクを使用して新しいマスク配列CPrasterを作成しています。

変更点のマスクの作成
- CL * CP_significance < cp_thres によって、信頼レベルと有意性の積が指定した閾値 cp_thres 未満の位置を見つけます。
- np.ma.mask_or を使用して、これらの位置と CL のマスクを組み合わせて新しいマスク cp_mask を作成します。
マスクのブロードキャスト
- cp_mask を Smasked の形状にブロードキャストします。これにより、Smasked と同じ形状のマスク cp_mask2 が作成されます。
新しいマスク配列の作成
- Smasked のデータを使用して、新しいマスク cp_mask2 を適用したマスク配列 CPraster を作成します。

To retrieve the dates of the change points we find the band indices in the time series along the time axis where the the maximum of the cumulative sums was located.
Numpy offers the "argmax" function for this purpose.

変化点の日付を取得するために、累積和の最大値が位置する時間軸に沿った時系列のバンドインデックスを見つける。
Numpyはこの目的のために "argmax "関数を提供している。

CP_index= np.ma.argmax(CPraster,axis=0)
change_indices = list(np.unique(CP_index))
change_indices.remove(0)
print(change_indices)
# Look up the dates from the indices to get the change dates
alldates=tindex[tindex>'2015-10-31']
change_dates=[str(alldates[x+1].date()) for x in change_indices]
print(change_dates)

変更点のインデックスとそれに対応する日付を求めています。

最も大きな値を持つインデックスを取得
- CPraster 内の各列（各ピクセル位置）で最も大きな値を持つインデックスを取得します。np.ma.argmax はマスクされた配列に対しても動作し、最も大きな値の位置を返します。
ユニークな変更インデックスのリストを作成し、0を削除
- np.unique を使用して、CP_index のユニークなインデックスを取得します。
- 変更インデックス 0 は無視されるため、それをリストから削除します。
インデックスから対応する日付を取得
- tindex から、2015年10月31日より後の日付を含む部分配列 alldates を取得します。
- 各変更インデックス x に対応する日付を alldates から取得し、文字列形式の日付のリスト change_dates を作成します。

Lastly, we visualize the change dates by showing the CP_index raster and label the change dates.

最後に、CP_indexラスタを表示し、変更日にラベルを付けることで、変更日を可視化する。

ticks=change_indices
ticklabels=change_dates

cmap=plt.cm.get_cmap('magma',ticks[-1])
fig, ax = plt.subplots(figsize=(8,8))
cax = ax.imshow(CP_index,interpolation='nearest',cmap=cmap)
# fig.subplots_adjust(right=0.8)
# cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
# fig.colorbar(p,cax=cbar_ax)

ax.set_title('Dates of Change')
# cbar = fig.colorbar(cax,ticks=ticks)
cbar=fig.colorbar(cax,ticks=ticks,orientation='horizontal')
_=cbar.ax.set_xticklabels(ticklabels,size=10,rotation=45,ha='right')

Secondary Change Points

After detection of a change point in the time series we can split the series in before and after change subsets. For forest degradation or deforestation detection for example this could apply when over the course of a multi-year time series selective logging precedes a clearing event or conversion of a logged plot to agriculture or regrowth, which show typically different time series profiles of radar backscatter. The approach to detect secondary change points would be to repeat analysis of the time series split into before and after change point detection.

時系列の変化点を検出した後は、変化前と変化後で時系列データを分割することができる。例えば、森林劣化や森林減少を検出する場合、数年にわたる時系列の中で、選択的伐採が皆伐に先行した場合や、伐採された区画が農業や再生林に転換された場合に適用できる。二次的な変化点を検出するアプローチは、変化点検出の前後に分割した時系列の分析を繰り返すことである。

conclusion

Pandas and numpy are powerful open source scripting tools to implement change point detection on large data stacks. For image based analysis numpy offers more efficient implementations compared to pandas, whereas pandas is more powerful in date time processing, e.g. time-weighted interpolation.

Pandas と numpy は、大規模なデータスタックの変化点検出を実装するための強力なオープンソーススクリプトツールです。画像ベースの分析では、numpyはpandasに比べてより効率的な実装を提供し、一方、pandasは、時間加重補間などの日付時間処理においてより強力です。

まとめ

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