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

はじめに

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

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

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 3 - Change Detection with Time Series Metrics and Log Ratio Method

Josef Kellndorfer, Ph.D., President and Senior Scientist, Earth Big Data, LLC
Revision date: January 2018
In this chapter we introduce three methods for change detection based on
Time series metrics 95 𝑡ℎ and 5 𝑡ℎ percentile difference thresholding
Time series coefficient of variation thresholding
Log Ratio from two image pairs

このPartでは、以下の3つの手法に基づく変化検出を紹介します。
・時系列データの95パーセンタイルと5パーセンタイルの差を用いたしきい値法
・時系列データの変動係数を用いたしきい値法
・2組の画像ペアを使った対数比法

Import Python modules

import os,sys,gdal
%matplotlib inline
import matplotlib.pylab as plt
import matplotlib.patches as patches  # Needed to draw rectangles
from skimage import exposure # to enhance image display
import numpy as np
import pandas as pd

まずはライブラリをインポートしています。

Select the Project data set and time series data

Louisiana Timber Management Site

# SENTINEL-1 TIME SERIES STACK VV from LOUISIANA FOREST MANAGEMENT SITE
datapath='/dev/shm/projects/c303nisar/louisiana/15SWRsS1/15SWRsS1-EBD/'
imagefile='15SWRsS1_A_vv_0063_A_mtfil.vrt'
datefile='15SWRsS1_A_vv_0063_A_mtfil.dates'

West Africa - Biomass Site

datapath='/Users/josefk/servir_data/wa/BIOsS1'
datefile='S32631X398020Y1315440sS1_A_vv_0001_mtfil.dates'
imagefile='S32631X398020Y1315440sS1_A_vv_0001_mtfil.vrt'
imagefile_cross='S32631X398020Y1315440sS1_A_vh_0001_mtfil.vrt'

os.chdir(datapath)

os.chdir(path)は、Pythonのosモジュールの関数で、現在の作業ディレクトリを指定されたパスに変更します。
datapathは、先ほど指定されたデータが格納されているディレクトリのパスです。（'/dev/shm/projects/c303nisar/louisiana/15SWRsS1/15SWRsS1-EBD/'）

We are defining two helper functions for this task
CreateGeoTiff() to write out images
timeseries_metrics() to compute various metrics from a time series data stack

今回、補助関数として以下の２つの関数を定義します。
・CreateGeoTiff()
目的：画像データを書き出すための関数
役割：地理空間データをGeoTIFF形式で保存します。GeoTIFFは、地理空間情報を持つTIFF形式のファイルで、GIS（地理情報システム）で広く使用されます。この関数は、処理した画像データや解析結果をファイルとして保存する際に使用されます。

・timeseries_metrics()
目的：時系列データスタックからさまざまなメトリクスを計算する関数
役割：与えられた時系列データから特定の統計情報やメトリクスを計算します。これは、変化検出やパターン認識のための重要なステップです。例えば、平均値、中央値、変動係数などの統計量を計算し、時系列データの特性を理解するために使用されます。

def CreateGeoTiff(Name, Array, DataType, NDV,bandnames=None,ref_image=None, 
                  GeoT=None, Projection=None):
    # If it's a 2D image we fake a third dimension:
    if len(Array.shape)==2:
        Array=np.array([Array])
    if ref_image==None and (GeoT==None or Projection==None):
        raise RuntimeWarning('ref_image or settings required.')
    if bandnames != None:
        if len(bandnames) != Array.shape[0]:
            raise RuntimeError('Need {} bandnames. {} given'
                               .format(Array.shape[0],len(bandnames)))
    else:
        bandnames=['Band {}'.format(i+1) for i in range(Array.shape[0])]
    if ref_image!= None:
        refimg=gdal.Open(ref_image)
        GeoT=refimg.GetGeoTransform()
        Projection=refimg.GetProjection()
    driver= gdal.GetDriverByName('GTIFF')
    Array[np.isnan(Array)] = NDV
    DataSet = driver.Create(Name, 
            Array.shape[2], Array.shape[1], Array.shape[0], DataType)
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection)
    for i, image in enumerate(Array, 1):
        DataSet.GetRasterBand(i).WriteArray( image )
        DataSet.GetRasterBand(i).SetNoDataValue(NDV)
        DataSet.SetDescription(bandnames[i-1])
    DataSet.FlushCache()
    return Name

この関数「CreateGeoTiff」は、与えられたデータをGeoTIFF形式で保存するためのものです。GeoTIFFは、地理空間情報を含むTIFFファイルの形式で、GIS（地理情報システム）で広く使われます。
以下に各プログラムの解説を記載します。
「CreateGeoTiff」のパラメータは、出力ファイル名(Name)、データ配列(Array)、データ型(DataType)、NoData値(NDV)、バンド名(bandnames)、参照画像(ref_image)、ジオトランスフォーム(GeoT)、投影情報(Projection)です。
各々のパラメータ配下を示しています。

• Name：出力ファイルの名前
• Array：書き出す画像データを含むNumPy配列であり、配列の形状は(バンド数, 高さ, 幅)の3次元配列
• DataType：GDALでサポートされるデータ型を指定（例：gdal.GDT_Float32等)
• NDV：データ内の欠損値を示す値を指定(NaNとして扱われる要素に対して設定)
• bandnames：各バンドの名前を指定するオプションのリスト
• ref_image：GeoTIFFファイルからジオトランスフォームおよび投影情報を取得するために使用するオプションの画像ファイル
  GeoT：画像の空間的な配置を定義するための6つの要素からなるタプル（例：(x_origin, pixel_width, 0, y_origin, 0, pixel_height)）
  Projection：画像の地理座標系を定義するための投影情報（WKT形式の文字列）

# If it's a 2D image we fake a third dimension:
    if len(Array.shape)==2:
        Array=np.array([Array])

上記のコードは入力データが2次元配列の場合、3次元配列に変換しています。これは、GDALが3次元配列（バンド x 高さ x 幅）としてデータを扱うためです。

    if ref_image==None and (GeoT==None or Projection==None):
        raise RuntimeWarning('ref_image or settings required.')

上記のコードは参照画像(ref_image)が指定されていない、かつジオトランスフォーム(GeoT)と投影情報(Projection) の必須チェックです。どちらも指定されていない場合、警告を発生させます。

    if bandnames != None:
        if len(bandnames) != Array.shape[0]:
            raise RuntimeError('Need {} bandnames. {} given'
                               .format(Array.shape[0],len(bandnames)))

上記のコードはバンド名が指定されている場合、バンド名の数がデータのバンド数と一致するかを確認しています。一致しない場合、エラーを発生させます。

    else:
        bandnames=['Band {}'.format(i+1) for i in range(Array.shape[0])]

上記のコードはバンド名が指定されていない場合、デフォルトで「Band 1」、「Band 2」などの名前を設定しています。

    if ref_image!= None:
        refimg=gdal.Open(ref_image)
        GeoT=refimg.GetGeoTransform()
        Projection=refimg.GetProjection()

上記のコードは参照画像が指定されている場合、その画像からジオトランスフォームと投影情報を取得しています。

    driver= gdal.GetDriverByName('GTIFF')

上記のコードはGeoTIFF形式のドライバを取得します。これを使って新しいGeoTIFFファイルを作成します。

    Array[np.isnan(Array)] = NDV

上記のコードはデータ配列内のNaN（Not a Number）値をNoData値(NDV)に置き換えています。

    DataSet = driver.Create(Name, 
            Array.shape[2], Array.shape[1], Array.shape[0], DataType)

上記のコードは指定された名前、データの幅、高さ、バンド数、データ型で新しいGeoTIFFファイルを作成しています。

    DataSet.SetGeoTransform(GeoT)

上記のコードはジオトランスフォーム情報を設定しています。これにより、画像の地理的な位置が定義されます。

    DataSet.SetProjection(Projection)

上記のコードは投影情報を設定しています。これにより、画像がどの地理座標系に属しているかが定義されます。

    for i, image in enumerate(Array, 1):
        DataSet.GetRasterBand(i).WriteArray( image )
        DataSet.GetRasterBand(i).SetNoDataValue(NDV)
        DataSet.SetDescription(bandnames[i-1])

上記のコードは各バンドに対して、データを書き込み、NoData値を設定し、バンド名を説明として設定しています。

    DataSet.FlushCache()

上記のコードはキャッシュされたデータをディスクに書き出しています。

    return Name

最終的に作成されたGeoTIFFファイルの名前を返しています。

次に「timeseries_metrics」関数の定義です。

def timeseries_metrics(raster,ndv=0): 
    # Make us of numpy nan functions
    # Check if type is a float array
    if not raster.dtype.name.find('float')>-1:
        raster=raster.astype(np.float32)
    # Set ndv to nan
    if ndv != np.nan:
        raster[np.equal(raster,ndv)]=np.nan
    # Build dictionary of the metrics
    tsmetrics={}
    rperc = np.nanpercentile(raster,[5,50,95],axis=0)
    tsmetrics['mean']=np.nanmean(raster,axis=0)
    tsmetrics['max']=np.nanmax(raster,axis=0)
    tsmetrics['min']=np.nanmin(raster,axis=0)
    tsmetrics['range']=tsmetrics['max']-tsmetrics['min']
    tsmetrics['median']=rperc[1]
    tsmetrics['p5']=rperc[0]
    tsmetrics['p95']=rperc[2]
    tsmetrics['prange']=rperc[2]-rperc[0]
    tsmetrics['var']=np.nanvar(raster,axis=0)
    tsmetrics['cov']=tsmetrics['var']/tsmetrics['mean']
    return tsmetrics

この関数は時系列ラスターデータから色々な統計メトリクスを計算しています。
以下にプログラムの解説を記載します。

def timeseries_metrics(raster, ndv=0):
    # Make use of numpy nan functions

本関数はラスターデータ(raster)とNoData値(ndv、デフォルトは0)を入力として受け取ります。

    # Check if type is a float array
    if not raster.dtype.name.find('float') > -1:
        raster = raster.astype(np.float32)

入力ラスターデータが浮動小数点型であるかを確認します。もしそうでなければ、float32型に変換します。

    # Set ndv to nan
    if ndv != np.nan:
        raster[np.equal(raster, ndv)] = np.nan

NoData値（ndv）が np.nan でない場合、ラスターデータ内のNoData値を np.nan に置き換えます。これにより、NaN関数（例えば、np.nanmean）を使用して計算が可能になります。

    # Build dictionary of the metrics
    tsmetrics = {}

計算結果を格納するための空の辞書「tsmetrics」を作成します。

    rperc = np.nanpercentile(raster, [5, 50, 95], axis=0)

ラスターデータの5パーセンタイル、50パーセンタイル（中央値）、95パーセンタイルを計算し、「rperc」に格納します。
パーセンタイルとは、データセットを百分位に分けた位置を示す値のことです。パーセンタイルは、データの分布を理解するために使用され、特定のパーセンタイル値はデータの特定の割合以下の値を示します。
以下に例示を記載します。
5パーセンタイル (p5)：データセットの中で、最小から数えて5％の位置にある値を示します。これは、データの下位5％に含まれる値です。
50パーセンタイル (p50)：データセットの中で、中央の位置（中央値）にある値を示します。
95パーセンタイル (p95)：データセットの中で、上位から数えて5％の位置にある値を示します。これは、データの95％がこの値以下にあることを意味します。

    tsmetrics['mean'] = np.nanmean(raster, axis=0)
    tsmetrics['max'] = np.nanmax(raster, axis=0)
    tsmetrics['min'] = np.nanmin(raster, axis=0)
    tsmetrics['range'] = tsmetrics['max'] - tsmetrics['min']
    tsmetrics['median'] = rperc[1]
    tsmetrics['p5'] = rperc[0]
    tsmetrics['p95'] = rperc[2]
    tsmetrics['prange'] = rperc[2] - rperc[0]
    tsmetrics['var'] = np.nanvar(raster, axis=0)
    tsmetrics['cov'] = tsmetrics['var'] / tsmetrics['mean']

上から順に、
平均値の計算：NaNを無視して、ラスターデータの平均値を計算し、「mean」として辞書に格納します。
最大値の計算：NaNを無視して、ラスターデータの最大値を計算し、「max」として辞書に格納します。
最小値の計算：NaNを無視して、ラスターデータの最小値を計算し、「min」として辞書に格納します。
範囲の計算：最大値と最小値の差を計算し、「range」として辞書に格納します。
中央値の設定：計算済みの50パーセンタイル値を「median」として辞書に格納します。
5パーセンタイルの設定：計算済みの5パーセンタイル値を「p5」として辞書に格納します。
95パーセンタイルの設定：計算済みの95パーセンタイル値を「p95」として辞書に格納します。
パーセンタイル範囲の計算：95パーセンタイル値と5パーセンタイル値の差を計算し、「prange」として辞書に格納します。
分散の計算：NaNを無視して、ラスターデータの分散を計算し、「var」として辞書に格納します。
変動係数の計算：分散を平均値で割った値を計算し、「cov」として辞書に格納します。

    return tsmetrics

最後に計算された結果を含む辞書「tsmetrics」を返します。

Set the Dates

# Get the date indices via pandas
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()

Bands and dates for S32631X398020Y1315440sS1_A_vv_0001_mtfil.vrt
1 2015-03-22 2 2015-04-03 3 2015-04-15 4 2015-05-09 5 2015-05-21
6 2015-06-02 7 2015-06-14 8 2015-06-26 9 2015-07-08 10 2015-07-20
11 2015-08-01 12 2015-08-13 13 2015-08-25 14 2015-09-06 15 2015-09-18
16 2015-09-30 17 2015-10-12 18 2015-10-24 19 2015-11-17 20 2015-11-29
21 2015-12-11 22 2015-12-23 23 2016-01-04 24 2016-01-28 25 2016-02-09
26 2016-03-04 27 2016-03-16 28 2016-03-28 29 2016-04-09 30 2016-04-21
31 2016-05-03 32 2016-05-15 33 2016-05-27 34 2016-06-08 35 2016-07-02
36 2016-07-14 37 2016-07-26 38 2016-08-07 39 2016-08-19 40 2016-08-31
41 2016-09-12 42 2016-09-24 43 2016-10-06 44 2016-10-18 45 2016-10-30
46 2016-11-11 47 2016-11-23 48 2016-12-05 49 2016-12-17 50 2016-12-29
51 2017-01-10 52 2017-01-22 53 2017-02-03 54 2017-02-15 55 2017-02-27
56 2017-03-11 57 2017-03-23 58 2017-04-04 59 2017-04-16 60 2017-04-28
61 2017-05-10 62 2017-05-22 63 2017-06-03 64 2017-06-15 65 2017-06-27
66 2017-07-09 67 2017-07-21 68 2017-08-02 69 2017-08-14 70 2017-08-26
71 2017-09-07 72 2017-09-19 73 2017-10-13 74 2017-10-25 75 2017-11-06
76 2017-11-18 77 2017-11-30

fi
このコードは、指定された日付ファイルから日付データを読み込み、各バンドと対応する日付を表示するためのものです。

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

この部分はPart2でも解説しているので割愛します。

for i in tindex:
    print("{:4d} {}".format(j, i.date()),end=' ')
    j+=1
    if j%5==1: print()

このコードは日付があるごとにバンド番号「j」と対応する日付「i.date()」をフォーマットして表示します。「{:4d}」はバンド番号を4桁で表示するためのフォーマット指定子です。バンド番号が5の倍数のたびに改行を挿入します。単純に出力整形をしているだけです。

Explore the Images

Below are a couple of plots showing the data set

以下はデータセットを示すいくつかのプロットです。

Open the image and get dimensions (bands,lines,pixels):

画像を開いて寸法(バンド、ライン、ピクセル)を取得します。

img=gdal.Open(imagefile)
img.RasterCount,img.RasterYSize,img.RasterXSize

(77, 3776, 4243)

GDAL（Geospatial Data Abstraction Library）を使って画像ファイルを開き、その画像に関する情報を取得するものです。具体的には、画像のバンド数、画像の高さ（ピクセル数）、画像の幅（ピクセル数）を取得します。

For a managable size for this notebook we choose a 1000x1000 pixel subset to read the entire data stack. We also convert the amplitude data to power data right away and will perform the rest of the calculations on the power data to be mathmatically correct. NOTE: Choose a different xsize/ysize in the subset if you need to.

以下のノートブックのサイズを管理しやすくするために、1000x1000ピクセルのサブセットを選択して、データスタック全体を読み込みます。また、振幅データをすぐにパワーデータに変換し、以降の計算はすべてパワーデータに対して行います。これにより、数学的に正確な処理が行われます。注意：必要に応じて、サブセットのサイズ（xsize/ysize）を変更してください。
要は大量データを一気に扱うと固まったり、予期せぬエラーを起こす可能性があるので1000×1000に一旦サイズを指定しますということです。

subset=(1500,0,500,500)   # (xoff,yoff,xsize,ysize)
bandnbr=1

rasterDN=img.GetRasterBand(bandnbr).ReadAsArray()
fig, ax = plt.subplots(figsize=(8,8))
ax.set_title('Sentinel-1 C-VV, NIGER!!!!!, {}'
             .format(tindex[bandnbr-1].date()))
ax.imshow(rasterDN,cmap='gray',vmin=2000,vmax=8000)
ax.grid(color='blue')
ax.set_xlabel('Pixels')
ax.set_ylabel('Lines')
# plot the subset as rectangle
if subset != None:
    _=ax.add_patch(patches.Rectangle((subset[0],subset[1]),
                                     subset[2],subset[3],
                                     fill=False,edgecolor='red',
                                     linewidth=3))

上記のコードはGDALとMatplotlibを使用して、特定のバンドの衛星画像データを読み込み、サブセットをハイライトして表示します。以下に、各行の詳細な説明を示します。

subset = (1500, 0, 500, 500)   # (xoff, yoff, xsize, ysize)
bandnbr = 1

表示する画像のサブセットを指定します。「subsetは(xoff, yoff, xsize, ysize)」の形式で、画像の左上角からのオフセット（xoff, yoff）、およびサブセットの幅（xsize）と高さ（ysize）を指定します。
また、表示するバンド番号を指定します。「bandnbr」はバンド番号を表し、今回は1番目のバンドを使用しています。

rasterDN = img.GetRasterBand(bandnbr).ReadAsArray()

指定したバンド番号の画像データを読み込み、NumPy配列 rasterDN に格納します。

fig, ax = plt.subplots(figsize=(8, 8))

Matplotlibを使用して、図(fig)と軸(ax)を作成します。図のサイズは8x8インチです。

ax.set_title('Sentinel-1 C-VV, NIGER!!!!!, {}'.format(tindex[bandnbr-1].date()))

画像のタイトルを設定します。ここでは、表示するバンドの取得日付を含むタイトルを設定します。

ax.imshow(rasterDN, cmap='gray', vmin=2000, vmax=8000)

読み込んだ画像データをグレースケールで表示します。表示範囲をvmin=2000からvmax=8000に設定しています。

ax.grid(color='blue')
ax.set_xlabel('Pixels')
ax.set_ylabel('Lines')

青色のグリッド線を表示し、X軸に「Pixels」、Y軸に「Lines」というラベルを設定します。

# plot the subset as rectangle
if subset != None:
    _ = ax.add_patch(patches.Rectangle((subset[0], subset[1]), subset[2], subset[3], fill=False, edgecolor='red', linewidth=3))

「subset」が「None」でない場合、サブセット領域を赤い長方形でハイライトします。patches.Rectangle を使用して、指定されたサブセットの位置とサイズで長方形を描画し、塗りつぶしなし(fill=False)、赤色の縁(edgecolor='red')、線幅3ピクセル(linewidth=3)で表示します。

rasterDN=img.ReadAsArray(*subset)
mask=rasterDN==0
CF=np.power(10.,-8.3)
rasterPwr=np.ma.array(np.power(rasterDN,2.)*CF,mask=mask,dtype=np.float32)
# Code below is an example to generate an 8bit scaled dB image
# rasterDB=(10.*np.ma.log10(rasterPwr)+31)/0.15   
# rasterDB[rasterDB<1.]=1.
# rasterDB[rasterDB>255.]=255.
# rasterDB=rasterDB.astype(np.uint8)
# rasterDB=rasterDB.filled(0)

上記コードはサブセット領域から画像データを読み込み、振幅データをパワーデータに変換しています。
コメントアウトされているコードは振幅データをデシベルスケールの8ビット画像に変換する例を示しています。
以下に、各行の詳細な説明を示します。

rasterDN = img.ReadAsArray(*subset)

GDALを使用して、指定されたサブセット領域 (subset) のデータを読み込みます。rasterDN には読み込まれたデータがNumPy配列として格納されます。

mask = rasterDN == 0

データ内の値が0である位置をマスクとして定義します。これは、パワーデータ計算時に0の値を無視するために使用されます。

CF = np.power(10., -8.3)

キャリブレーションファクターを定義します。この例では、振幅データをパワーデータに変換するため「10^-8.3」を使用します。

rasterPwr = np.ma.array(np.power(rasterDN, 2.) * CF, mask=mask, dtype=np.float32)

振幅データ(rasterDN)の平方を取ってパワーデータに変換し、キャリブレーションファクターを掛けます。マスクを適用して、値が0の位置を無視し、結果をfloat32型で格納します。

# Code below is an example to generate an 8bit scaled dB image
# rasterDB = (10. * np.ma.log10(rasterPwr) + 31) / 0.15

パワーデータをデシベルスケールに変換します。ここでは、デシベル値を計算し、その後8ビットスケールに変換します。

# rasterDB[rasterDB < 1.] = 1.
# rasterDB[rasterDB > 255.] = 255.

デシベルスケールの値が1未満の場合、1に設定し、デシベルスケールの値が255を超える場合、255に設定します。

# rasterDB = rasterDB.astype(np.uint8)

データを8ビット符号なし整数(uint8)に変換します。

# rasterDB = rasterDB.filled(0)

マスクされた値を0に置き換えます。これにより、データ内の欠損値が0で埋められます。

We make an RGB stack to display the first, center, and last time step as a multi-temporal color composite. The np.dstack results in an array of the form [lines,pixels,bands], which is the format we need for RGB display with matplotlib's imshow() function. Note that numpy array indexing starts with 0, so band 1 is raster[0].

最初の時点、中間の時点、最後の時点を使って、カラー画像として表示するためにRGBスタックを作ります。np.dstack を使うと、データが [行, 列, バンド] の形式になります。これはMatplotlibの imshow() 関数でRGB表示するために必要な形式です。注意点として、NumPyの配列は0から始まるので、バンド1は raster[0] になります。

rgb_bands=(1,int(img.RasterCount/2),img.RasterCount)  # first, center, last band
rgb_bands=(1,10,40) 
rgb_bands=(18,45,74)
rgb_idx=np.array(rgb_bands)-1  # get array index from bands by subtracting 1
rgb=np.dstack((rasterPwr[rgb_idx[0]],rasterPwr[rgb_idx[1]],rasterPwr[rgb_idx[2]]))
rgb_dates=(tindex[rgb_idx[0]].date(),
           tindex[rgb_idx[1]].date(),tindex[rgb_idx[2]].date())

このコードは、指定されたバンドを使ってRGB画像を作成し、各バンドに対応する日付を取得しています。以下に、実施していることをリストします。

• RGBバンドの選択
  最初の行では、最初のバンド、中間のバンド、最後のバンドを選択しています。
  次の行では、具体的なバンド番号（1, 10, 40）を指定しています。
  さらにその次の行では、他の具体的なバンド番号（18, 45, 74）を指定しています。
• インデックスの調整
  指定したバンド番号を配列インデックスに変換するため、各バンド番号から1を引いてインデックスを取得しています。
• RGBスタックの作成
  選択したバンドのパワーデータを使って、RGB画像データを作成しています。np.dstack を使用して3つのバンドをスタックし、RGB形式の配列を生成します。
• 日付の取得
  各選択したバンドに対応する日付を取得し、rgb_dates に格納しています。

We are also interested in displaying the image enhanced with histogram equalization. We can use the function exposure.equalize_hist() from the skimage.exposure module

画像をヒストグラム平坦化で見やすく表示したいと考えています。これには、skimage.exposure モジュールの exposure.equalize_hist() 関数を使います。
ヒストグラム平坦化（ヒストグラム均等化）とは、画像の輝度値の分布を均等にすることで、暗い部分を明るくし、明るい部分を暗くするなどして、画像全体のコントラストを高めます。

rgb_stretched=rgb.copy()
# For each band we apply the strech
for i in range(rgb_stretched.shape[2]):
    rgb_stretched[:,:,i] = exposure.\
    equalize_hist(rgb_stretched[:,:,i].data,
    mask=~np.equal(rgb_stretched[:,:,i].data,0.))

このコードは、RGB画像の各バンドに対してヒストグラム平坦化を適用して、コントラストを強調する処理を行っています。
以下に、実施していることをリストします。

• RGB画像のコピーを作成
  元のRGB画像をコピーして rgb_stretched に格納します。このコピーに対してヒストグラム平坦化を適用します。
• 各バンドにヒストグラム平坦化を適用
  • ループの開始
    RGB画像の各バンドに対して処理を行うために、ループを開始します。
    rgb_stretched.shape[2] はバンド数を示します。
  • ヒストグラム平坦化の適用
    「exposure.equalize_hist()」関数を使用して、各バンドにヒストグラム平坦化を適用します。
    • rgb_stretched[:, :, i].data
      現在のバンドのデータを取得します。
    • mask=~np.equal(rgb_stretched[:, :, i].data, 0.)
      マスクを作成します。ここでは、値が0でないピクセルに対してヒストグラム平坦化を適用します。

Now let's display the the unstrechted and histogram equalized images side by side.

では、ストレッチしていない画像とヒストグラム平坦化した画像を並べて表示してみましょう。

fig,ax = plt.subplots(1,2,figsize=(16,8))
fig.suptitle('Multi-temporal Sentinel-1 backscatter image R:{} G:{} B:{}'
             .format(rgb_dates[0],rgb_dates[1],rgb_dates[2]))
plt.axis('off')
ax[0].imshow(rgb)
ax[0].set_title('Unstreched')
ax[0].axis('off')
ax[1].imshow(rgb_stretched)
ax[1].set_title('Histogram Equalized')
_=ax[1].axis('off')

上記コードは図の表示コードであり、今までのコードと変わりないので簡単にやっていることをリストします。

• 2つのサブプロットを持つ図を作成
• 図全体のタイトルに各RGBバンドの日付を表示
• 左側のサブプロットにストレッチしていない画像を表示
• 右側のサブプロットにヒストグラム平坦化した画像を表示
• 両方の画像の軸を非表示にして、画像だけを見やすく表示

Computation and Visualization of the Time Series Metrics

For the entire time series, we will compute some metrics that will aid us in change detection. For each pixel in the stack we compute:
Mean
Median
Maximum
Minimum
Range (Maximum - Minimum)
5th Percentile
95th Percentile
PRange (95th - 5th Percentile)
Variance
Coefficient of Variation (Variance/Mean)

全ての時系列データに対して、変化検出に役立ついくつかのメトリクスを計算します。各ピクセルについて以下のメトリクスを計算します。

• 平均値
• 中央値
• 最大値
• 最小値
• 範囲（最大値 - 最小値）
• 5パーセンタイル
• 95パーセンタイル
• パーセンタイル範囲（95パーセンタイル - 5パーセンタイル）
• 分散
• 変動係数（分散 / 平均値）

metrics=timeseries_metrics(rasterPwr.filled(np.nan),ndv=np.nan)

このコードは、timeseries_metrics 関数を使って、時系列のパワーデータからメトリクスを計算しています。
パラメータについて補足解説します。
「rasterPwr.filled(np.nan)」はマスクされている値（欠損値）を「np.nan」で埋めます。これにより、欠損値がNaNとして扱われます。
「ndv=np.nan」は関数に渡すNoData値として「np.nan」を指定しています。

metrics.keys()

dict_keys(['mean', 'max', 'min', 'range', 'median', 'p5', 'p95', 'prange', 'va
r', 'cov'])

「metrics」に含まれる全てのメトリクスの名前（キー）をリストとして取得できます。

Let's look at the histograms for the time series variance and coeficient of variation to aid displaying those images:

時系列の分散と変動係数の画像を表示するために、それぞれのヒストグラムを確認してみましょう。

fig, ax= plt.subplots(1,2,figsize=(16,4))
ax[0].hist(metrics['var'].flatten(),bins=100)
ax[1].hist(metrics['cov'].flatten(),bins=100)
_=ax[0].set_title('Variance')
_=ax[1].set_title('Coefficient of Variation')

このコードは、以下のことを行っています。

• 2つのサブプロットを持つ図を作成
• 左のサブプロットに分散のヒストグラムを表示
• 右のサブプロットに変動係数のヒストグラムを表示
• 各サブプロットに適切なタイトルを設定

We use thresholds determined from those histograms to set the scaling in the time series visualiztion. For the backscatter metrics we choose a typeical range appropriate for this ecosystem and radar sensor. A typical range is -30 dB (0.0001) to -5.2 dB (0.3).

これらのヒストグラムから適切な閾値を決めて、時系列画像の表示範囲を設定します。散乱強度については、この環境とレーダーセンサーに合った範囲を選びます。一般的には「-30dB(0.0001)」から「-5.2dB(0.3)」の範囲が適しています。

# List the metrics keys you want to plot
metric_keys=['mean', 'median', 'max', 'min', 
             'p95', 'p5','range', 'prange','var','cov']
fig= plt.figure(figsize=(16,40))
idx=1
for i in metric_keys:
    ax = fig.add_subplot(5,2,idx)
    if i=='var': vmin,vmax=(0.0,0.005)
    elif i == 'cov': vmin,vmax=(0.,0.04)
    else:
        vmin,vmax=(0.0001,0.3)
    ax.imshow(metrics[i],vmin=vmin,vmax=vmax,cmap='gray')
    ax.set_title(i.upper())
    ax.axis('off')
    idx+=1

このコードは、指定されたメトリクスをプロットして表示するためのものです。
以下に、実施していることをリスト形式で解説します。

• プロットするメトリクスのキーをリストとして定義
• 16x40インチの図を作成
• 各メトリクスについてループ処理を行い、サブプロットを追加
• 各メトリクスに適した表示範囲を設定
  • 分散(var)の場合(0.0, 0.005)
  • 変動係数(cov)の場合(0.0, 0.04)
  • それ以外のメトリクスは(0.0001, 0.3)
• メトリクスデータをグレースケールでプロット
• 各サブプロットのタイトルを設定
• サブプロットの軸を非表示
• 次のサブプロットのインデックスを更新し、上記処理を再起実行

Change detection with the Percentile Difference Threshold Method

In this method we find thresholds on the 95 𝑡ℎ and 5 𝑡ℎ percentile difference. The advantage to look at percentiles verus maximum minus minimum is that outliers and extremas in the time series are not influencing the result.
For our example, the historgram of the 95 𝑡ℎ and 5 𝑡ℎ percentile difference image looks like this:

この方法では、95パーセンタイルと5パーセンタイルの差分に基づいて閾値を見つけます。最大値と最小値の差を見るのに比べて、パーセンタイルを使う利点は、時系列データの外れ値や極端な値が結果に影響を与えないことです。
具体例として、95パーセンタイルと5パーセンタイルの差分画像のヒストグラムは次のようになります。

plt.hist(metrics['range'].flatten(),bins=100,range=(0,0.3))
_=plt.axvline(0.27,color='red')

このコードは、「metrics['range']」データのヒストグラムを表示し、特定の閾値を赤い縦線で示しています。

• metrics['range'].flatten()
  「metrics」の「range」キーに対応するデータ（範囲）の配列を平坦化（1次元配列に変換）します。
• bins=100
  ヒストグラムのビン（棒)の数を100に設定します。
• range=(0, 0.3)
  ヒストグラムの範囲を0から0.3に設定します。
  これにより、metrics['range']のデータ分布がヒストグラムとして表示されます。
  

Let's visualize the change pixels (cp) where the 95th - 5th percentile difference in the time series for each pixel (x,y) exceed a threshold 𝑡

cp_{x,y} = P_{x,y}^{95th} - P_{x,y}^{5th} > t

With 𝑡=0.15 the image looks like:

各ピクセル (x,y) における時系列データの95パーセンタイルと5パーセンタイルの差が閾値 𝑡 を超える変化ピクセル (cp) を視覚化しましょう。
閾値𝑡=0.15の場合、画像は次のようになります。

thres=0.25
plt.figure(figsize=(8,8))
mask=metrics['range']<thres # For display we prepare the inverse mask
maskpdiff=~mask # Store this for later output
plt.imshow(mask,cmap='gray')
plt.legend(['$p_{95} - p_5 > 0.15$'],loc='center right')
_=plt.title('Threshold Classifier on Percentile Difference ($P_{95} - P_5 > 0.15$)')

上記コードは指定された閾値に基づいて変化ピクセルを視覚化し、結果をプロットします。
以下に、実施していることをリストします。

• 変化を検出するための閾値(0.25)を設定
• 指定された閾値に基づいて、変化を検出するためのマスクを作成
• マスクを反転して逆マスクを作成し、後で使用するために保存
• 反転前のマスクをグレースケールで表示
• 凡例を設定
• プロットのタイトルを設定(0.15とされているが、プログラム例はthres=0.25のため0.25の表記が正しい。)

Change Detection with the Coefficient of Variation Method

We can set a threshold 𝑡 for the coefficient of variation image to classify change in the time series:

cp_{x,y} = \frac{\sigma_{x,y}^2}{\overline{X_{x,y}}} > t

Let's look at the histogram of the coefficient of variation:

時系列データの変化を分類するために、変動係数画像の閾値𝑡を設定できます。

plt.hist(metrics['cov'].flatten(),bins=100,range=(0,0.05))
_=plt.axvline(0.025,color='red')

上記コードは記載済みのコードと同じなので解説は割愛します。

With a threshold t=0.01 the change pixels would look like the following image:
閾値t=0.01を使用すると、変化ピクセルは次のような画像になります。

thres=0.025
mask=metrics['cov'] < thres
maskcv=~mask
plt.figure(figsize=(8,8))
plt.imshow(mask,cmap='gray')
_=plt.title('Threshold Classifier on Time Series Coefficient of Variation')

Change Detection with the Log Ratio Method

We compare two images from the same season in different years. First we look at global means of the backscatter images in the subset building a time series object of acquisition dates and global image means of backscatter.

異なる年の同じ季節に撮影された2つの画像を比較します。まず、サブセット内の散乱強度画像のグローバル平均値を取得し、取得日と散乱強度のグローバル平均値の時系列オブジェクトを作成します。

tsmean=10*np.log10(np.nanmean(rasterPwr.filled(np.nan),axis=(1,2)))

このコードは、各時点での散乱強度画像の平均値をデシベル（dB）スケールで計算し、時系列として取得します。
以下に、各パラメータの説明を示します。

• rasterPwr.filled(np.nan)
  「rasterPwr」はマスク付き配列で、マスクされている値を「np.nan（NaN）」で埋めます。これにより、欠損値がNaNとして扱われます。
• np.nanmean(rasterPwr.filled(np.nan), axis=(1, 2))
  各時点（各バンド）の散乱強度画像の平均値を計算します。
  axis=(1, 2)により、ピクセルレベルでの平均値を計算します。
  次元的に(バンド、画像縦方向、画像横方向)となっており、(1, 2)は各バンドに対する縦×横の画像全体を示しています。
• np.nanmean
  NaN値を無視して計算します。
• 10*np.log10
  計算された平均値をデシベル（dB）スケールに変換します。

We make a time series object to list the dates, mean backscatter in dB and band index number for the rasterPwr array:

取得日、散乱強度の平均値（dB単位）、および rasterPwr 配列のバンドインデックス番号をリストするために、時系列オブジェクトを作成します。

ts = pd.Series(tsmean,index=tindex)
for i in range(len(ts)):
    print(i,ts.index[i].date(),ts[i])

0 2015-03-22 -9.773781
1 2015-04-03 -9.814333
2 2015-04-15 -9.84827
3 2015-05-09 -10.075288
4 2015-05-21 -9.987606
5 2015-06-02 -9.835003
6 2015-06-14 -10.412914
7 2015-06-26 -10.64331
8 2015-07-08 -9.98234
9 2015-07-20 -9.159636
10 2015-08-01 -7.678219
11 2015-08-13 -8.60141
12 2015-08-25 -7.6070075
13 2015-09-06 -7.645421
14 2015-09-18 -6.655918
15 2015-09-30 -8.7717705
16 2015-10-12 -9.348694
17 2015-10-24 -9.547744
18 2015-11-17 -10.2138815
19 2015-11-29 -11.099142
20 2015-12-11 -11.029471
21 2015-12-23 -11.332901
22 2016-01-04 -11.346351
23 2016-01-28 -11.197915
24 2016-02-09 -11.145014
25 2016-03-04 -10.9366045
26 2016-03-16 -11.114582
27 2016-03-28 -11.000681
28 2016-04-09 -10.456753
29 2016-04-21 -11.031124
30 2016-05-03 -11.042203
31 2016-05-15 -11.248089
32 2016-05-27 -10.781347
33 2016-06-08 -10.7717905
34 2016-07-02 -10.622729
35 2016-07-14 -10.262638
36 2016-07-26 -9.969166
37 2016-08-07 -9.227007
38 2016-08-19 -8.372538
39 2016-08-31 -7.8771267
40 2016-09-12 -9.163029
41 2016-09-24 -9.04641
42 2016-10-06 -10.078144
43 2016-10-18 -10.534364
44 2016-10-30 -11.044583
45 2016-11-11 -11.120414
46 2016-11-23 -11.056729
47 2016-12-05 -11.187023
48 2016-12-17 -11.514052
49 2016-12-29 -11.376835
50 2017-01-10 -11.243304
51 2017-01-22 -11.204616
52 2017-02-03 -11.176929
53 2017-02-15 -11.093778
54 2017-02-27 -11.04459
55 2017-03-11 -10.92975
56 2017-03-23 -10.895084
57 2017-04-04 -11.270055
58 2017-04-16 -11.106432
59 2017-04-28 -11.091718
60 2017-05-10 -11.196309
61 2017-05-22 -10.516581
62 2017-06-03 -11.056223
63 2017-06-15 -10.081734
64 2017-06-27 -10.121447
65 2017-07-09 -10.04177
66 2017-07-21 -9.837778
67 2017-08-02 -9.248106
68 2017-08-14 -8.226735
69 2017-08-26 -8.096363
70 2017-09-07 -9.48612
71 2017-09-19 -9.06922
72 2017-10-13 -10.274279
73 2017-10-25 -10.869721
74 2017-11-06 -11.154692
75 2017-11-18 -11.122526
76 2017-11-30 -11.273689

このコードは、時系列データをpandasのSeriesオブジェクトとして作成し、各時点のインデックス（日付）と対応する散乱強度の平均値（dB）を表示します。以下に、各行の説明をリストします。

• 時系列オブジェクトの作成
  • tsmean
    以前に計算した散乱強度の平均値（dB）の配列です。
  • tindex
    各時点の日付を含むDatetimeIndexです。
  • pd.Series(tsmean, index=tindex)
    「tsmean」の値を「tindex」の日付に対応させたpandasのSeriesオブジェクト「ts」を作成します。これにより、日付と対応する平均値を持つ時系列データが得られます。
• 時系列データの表示
  tsの数だけループし、日付と散乱強度を表示します。

To compare two dates for change detection with the log ratio approach we pick two dates of relative low backscatter (dry conditions) and from similar times of the year. Two such candidate dates are

West Africa / Biomass Site example:
2015-11-29 -11.099142 dB (index 19)
2017-11-30 -11.273689 dB (index 76)

ログ比アプローチを用いて変化検出を行うために、散乱強度が比較的低い（乾燥条件）2つの日付を選び、さらに同じ時期に撮影されたものを選びます。
ログ比アプローチとは、2つの時点で取得された散乱強度データの比を対数を使って計算し、変化を検出する方法です。

# # Louisiana
# Xr=rasterPwr[7]   # Reference Image
# Xi=rasterPwr[37]  # New Image
# WA biomass
Xr=rasterPwr[19]   # Reference Image
Xi=rasterPwr[76]  # New Image

このコードは、特定の時点の散乱強度データを参照画像と新しい画像として設定しています。ルイジアナと西アフリカのバイオマスサイトの2つの例が示されていますが今回使用されているのは西アフリカのバイオマスになります。

The Log ratio between the images is

画像間の対数比率は次の式にて表せられます。

r  = log_{10}(\frac{X_i}{X_r})

r = np.log10(Xi/Xr)

To find a threshold for change, we can display the absolute ration image 𝑎𝑏𝑠(𝑟) and the historgram of 𝑟 . We adjust the scale factors for the display to enhance visualization of change areas with largest backscatter change over the time series. Brighter values show larger change.

変化を検出するための閾値を見つけるために、絶対比画像abs(r)と比率(r)のヒストグラムを表示します。
時系列データにおける散乱強度の大きな変化を視覚的に強調するために、表示のスケールを調整します。明るい値は、より大きな変化を示します。

# Display r
fig, ax = plt.subplots(2,1,figsize=(8,16))
ax[0].axis('off')
ax[0].imshow(np.abs(r),vmin=0,vmax=0.3,cmap='gray')
_=ax[1].hist(r.flatten(),bins=100,range=(-0.4,0.4))

上記コードは解説済みコードと同様なので説明は割愛します。

Let's define change pixels as those falling outside the range of three times the standard deviation of the ration image 𝜎𝑟 from the image mean 𝑟¯ :

{cp}_{x,y} = (r_{x,y} < \overline{r} - 3\sigma_r) \ \textrm{or} \ (r_{x,y} > \overline{r} + 3\sigma_r)

We are using the numpy masking to set the non-changing pixels inside the range:

変化ピクセルを、比率画像の平均𝑟‾から3倍の標準偏差3𝜎𝑟を超えて外れるピクセルとして定義しましょう。
この範囲内の変化しないピクセルを設定するために、NumPyのマスキングを使用します。

stddev=np.std(r)
thres=3*stddev
mask=np.logical_and(r>-1*thres,r<thres)
masklr=~mask

このコードは、比率画像 r に対して、3標準偏差を閾値として変化ピクセルをマスクするための準備を行います。
実施していることを簡単にリストすると以下になります。

• 比率画像rの標準偏差を計算
• 標準偏差の3倍を閾値として設定
• 比率画像の値が-3標準偏差から3標準偏差の範囲内にあるピクセルを示すマスクを作成
• そのマスクを反転して、変化ピクセルを示すマスクを作成

Let's display pixels that fall outside 3 times the standard deviation
3倍の標準偏差から外れたピクセルを表示してみましょう。

fig,ax = plt.subplots(figsize=(8,16))
ax.imshow(mask,cmap='gray')
ax.xaxis.set_ticks([])
ax.yaxis.set_ticks([])
_=ax.set_title('Log Ratio Classifier of the October 2016/2017 Log Ratio Images')

Write the images to an output file

Determine output geometry
First, we need to set the correct geotransformation and projection information. We retrieve the values from the input images and adjust by the subset:

出力のジオメトリを決定します。
まず、正しいジオトランスフォーメーションと投影情報を設定する必要があります。入力画像からこれらの値を取得し、サブセットに合わせて調整します。

proj=img.GetProjection()
geotrans=list(img.GetGeoTransform())

subset_xoff=geotrans[0]+subset[0]*geotrans[1]  
subset_yoff=geotrans[3]+subset[1]*geotrans[5] 
geotrans[0]=subset_xoff
geotrans[3]=subset_yoff
geotrans=tuple(geotrans)
geotrans

(428020.0, 20.0, 0.0, 1390960.0, 0.0, -20.0)

このコードは、入力画像のジオトランスフォーメーションと投影情報を取得し、サブセットに合わせて調整しています。以下に各行の詳細な説明を示します。

• 入力画像から投影情報を取得
• 入力画像からジオトランスフォーメーション情報を取得
• サブセットのオフセットを考慮して、新しいXオリジンとYオリジンを計算
  元のジオトランスフォーメーションのXオリジンに、サブセットのXオフセット（ピクセル単位）をピクセルサイズ（X方向）で掛け合わせた値を足して、新しいXオリジンsubset_xoffを計算します。Yも同様です。
• 新しいXオリジンとYオリジンをジオトランスフォーメーションに設定
• ジオトランスフォーメーションをタプルに変換

Time series metrics images

We use the root of the time series data stack name and append a ts_metrics.tif ending as filenames

時系列データスタックの名前に基づいて、ファイル名の末尾に "ts_metrics.tif" を追加した名前を使用します。

# Time Series Metrics as image image:
# We make a new subdirectory where we will store the images
dirname=imagefile.replace('.vrt','_tsmetrics2')
os.makedirs(dirname,exist_ok=True)
print(dirname)

S32631X398020Y1315440sS1_A_vv_0001_mtfil_tsmetrics2

このコードは、時系列メトリクスを画像として保存するために、新しいサブディレクトリを作成し、そのディレクトリ名を出力しています。実施内容をリストします。

• ディレクトリ名の生成
  • imagefile からファイル拡張子 .vrt を _tsmetrics2 に置き換えた名前を dirname に格納
    例えば、imagefile が example.vrt であれば、dirname は example_tsmetrics2 になります。
• サブディレクトリの作成
  • os.makedirs 関数を使って、dirname で指定された名前のサブディレクトリを作成
    exist_ok=True は、すでにディレクトリが存在する場合にエラーを発生させないようにするオプションです。

Output the individual metrics as GeoTIFF images:

Names=[] # List to keep track of all the names
for i in metrics:
    # Name, Array, DataType, NDV,bandnames=None,ref_image
    Name=os.path.join(dirname,imagefile.replace('.vrt','_'+i+'.tif'))
    CreateGeoTiff(Name,metrics[i],gdal.GDT_Float32,np.nan,[i],GeoT=geotrans,Projection=proj)
    Names.append(Name)

このコードは、すべての時系列メトリクスを含むGeoTIFFファイルを作成し、それらのファイル名をリストに保持します。以下に、各行の詳細な説明を示します。

• ファイル名のリストの初期化
  後で作成するすべてのGeoTIFFファイルの名前を格納するためのリスト Names を初期化します。
• メトリクスのループ
• ファイルの生成
  • 元の imagefile の拡張子「.vrt」を、現在のメトリクスキー「i」と「.tif」拡張子に置き換えた名前を生成
  • 生成されたファイル名は、前に作成したサブディレクトリ「dirname」に結合
    例えば、「imagefile」が「example.vrt」で、メトリクスキーが「mean」の場合、生成されるファイル名は「example_mean.tif」になります。
• GeoTIFFファイルの作成
• ファイル名の追加
  作成したGeoTIFFファイルの名前をリスト「Names」に追加します。

Build a Virtual Raster Table on the Metrics GeoTIFF images

To tie the images in to one new raster stack of time series metrics we build a virtual raster table with all the metrics.
Trick: Use ' '.join(Names) to build one long string of names separated by a space as input to gdalbuildvrt

時系列メトリクスの新しいラスタースタックを作成するために、すべてのメトリクスを含む仮想ラスターテーブルを構築します。
コツ：' '.join(Names)を使用して、名前をスペースで区切った長い文字列を作成し、gdalbuildvrtの入力として使用します。

cmd='gdalbuildvrt -separate -overwrite -vrtnodata nan '+\
   dirname+'.vrt '+' '.join(Names)
# print(cmd)
os.system(cmd)

0

このコードは、複数のGeoTIFFファイルを一つの仮想ラスタースタック（VRTファイル）に統合するためのコマンドを生成し、実行します。以下に、各行の詳細な説明を示します。

• 複数のGeoTIFFファイルを一つの仮想ラスタースタック（VRTファイル）に統合するためのgdalbuildvrt コマンドを生成
• 生成されたコマンドをシステムシェルで実行し、VRTファイルを作成

os.getcwd()

'/Users/rmuench/Downloads/wa/BIOsS1'

現在の作業ディレクトリを表示します。

print('Time Series Metrics VRT File:\n',dirname+'.vrt')

Time Series Metrics VRT File:
S32631X398020Y1315440sS1_A_vv_0001_mtfil_tsmetrics2.vrt

Change Images from the three methods

We are going to write one three-band GeoTIFF output file that stores the results from the three classifiers

3つの分類器の結果を保存するために、3バンドのGeoTIFF出力ファイルを作成します。

imagename=imagefile.replace('.vrt','_thresholds.tif')
bandnames=['Percentile','COV','Log Ratio']
Array=np.array([maskpdiff,maskcv,masklr])
CreateGeoTiff(imagename,Array,gdal.GDT_Byte,0,bandnames,GeoT=geotrans,Projection=proj)

このコードは、3つの分類器の結果を3バンドのGeoTIFFファイルとして保存します。リストすると以下になります。

• 出力ファイル名を生成
• 3つの分類器の名前をバンド名として定義
• 3つのマスクを含むデータ配列を作成
• CreateGeoTiff関数を使用して、3バンドのGeoTIFFファイルを作成し、3つの分類器の結果を保存

This Image can now be loaded into QGIS or other Tools and only the detected layers should show.

この画像はQGIS等のツールにロードできるようになり、検出されたレイヤーのみが表示されます。

Conclusion

Thresholds for the three methods are necessarily site dependent and need to be identified with calibration data or visual post-classification interpretation, and can subsequently be adjusted to maximize classification accuracy. Also, some methods will have advantages in different scenarios.
At the Earth Big Data SEPPO Processor we actually transform many of the time series metrics data types back to lower volume storage models, e.g. 16 bit scaled amplitudes. See the EBD Data Guide here(also available as pdf).

3つの方法の閾値は、場所によって異なるため、キャリブレーションデータや視覚的な後分類解釈を使って特定する必要があります。また、分類精度を最大化するために、閾値を調整することができます。さらに、特定のシナリオでは、各方法に異なる利点があります。

Earth Big DataのSEPPOプロセッサでは、実際に多くの時系列メトリクスデータ型を、16ビットのスケール振幅などの低容量ストレージモデルに変換しています。詳細は、EBDデータガイドをご参照ください（PDFでも利用可能）。

まとめ

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