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

はじめに

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

SERVIRとは以下のように紹介されています。

SERVIR is a NASA and USAID partnership that supports locally led efforts to strengthen climate resilience, food and water security, forest and carbon management, and air quality.

つまり、SERVIRはNASAとUSAIDのパートナーシップであり、気候変動に対する回復力、食糧と水の安全保障、森林と炭素の管理、大気の質を強化するための地元主導の取り組みを支援しているそうです。

本記事では「SAR Handbook」のChp3の練習問題について記載します。

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 1 - Getting to Know SAR Images and Forest Signatures

Josef Kellndorfer, Ph.D., President and Senior Scientist, Earth Big Data, LLC
Revision date: January 2018
This book chapter on SAR data analysis for forest applications focused on deforestation and forest degradation monitoring is implemented as an interactive notebook. The digital format (a jupyter notebook) of this chapter can readily be launched in any webbrowser for interactive data exploration with provided or new traning data. The notebook is comprised of text written in a combination of executable python code and markdown formatting including latex style mathmatical equations. With this approach, the trainees can readily expand, change and share the entire work with new data sets in new regions or newly available time series steps.
While we are only scratching the surface of available open source tools, the course will provide a broad overview on what modern tools can be employed for focused SAR data analysis, or remote sensing data analysis in general.

ここで紹介されている問題はJupyter Notebook形式で配布しているため任意のWebブラウザで実行できます。

Software Installation and Data Sets
Please refer to the documents INSTALLATION and DATA_HOWTO.
The time series data sets for this training course were pre-processed with the EARTH BIG DATA Software for Earth Big Data Processing, Prediction Modeling and Organization (SEPPO) using cloud-based processing on Amazon Web Services. SEPPO allows for the fully automated processing of large SAR (and other remote sensing) data sets to constuct time series data effectively. The data format guide EBD_README explains data structures and filenaming conventions for data sets produced by EARTH BIG DATA, LLC.

本章で扱うデータはAWS(Amazon Web Services)のサービスによって前処理されたデータを使用します。

Importing relevant python packages

First step in the time series analysis approach after obtaining the preprocessed data stacks is the import of necessary python packages. See the comments below as to what packages are needed and their functions. Note that all these packages should have been installed when the python anaconda environment was created.

前処理されたデータを取得したら、時系列分析する際の最初のステップは、必要なpythonパッケージのインポートです。必要なパッケージは、後述します。これらのパッケージはすべて、python anaconda環境が作成されたときにインストールされているはずです。

import pandas as pd
import gdal
import numpy as np
import time,os, glob
%matplotlib inline
import matplotlib.pylab as plt

pandas はデータ操作や分析のためのライブラリです。
gdal は地理空間データを操作するためのライブラリです。
numpy は数値計算を効率的に行うためのライブラリです。
time, os, glob はPythonの標準ライブラリで、ファイル操作や時間計測などに使います。
matplotlibは描画用のライブラリです。
%matplotlib inlineとすることによって、Jupyter上、つまりブラウザ上でも表示できるようにしています。

import gdalでそんなライブラリ無いとエラーになった人は以下のimportに直して実行してみてください。
from osgeo import gdal

Set Project Directory and Filenames

Edit and uncomment the respective cell entries below to activate the wanted project data directory. Take a look at the EBD Data Guide for an explanation of the filenaming conventions used for image and date files.
How to specify data directories:
Linux: /path/to/file
Windows: d:/path/to/file
d: is the drive letter # IMPORTANT: Always use '/' instead '' in Windows
NOTE: directories and filenames are specified in python as strings enclosed in single or double quotes: 'string' "string"

データ・ディレクトリの格納場所を指定するため、以下のようにパスを指定してください。
画像ファイルと日付ファイルのファイル名の規則については、EBDデータガイドを参照してください。
データディレクトリの指定方法の例示は以下です。
Linux：/パス/to/ファイル
Windows: d:/path/to/file
注意: ディレクトリとファイル名は、Pythonでは一重引用符(')または二重引用符(")で囲んでください。

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'

4種類分のデータファイルについて、ディレクトリとファイル名の指定の例示になります。
ちなみに筆者はよくGoogleDriveを利用していますが、Driveに格納した場合、GoogleColabからDriveを参照し、それぞれのファイルを右クリックするとパスがコピーできます。
また、データは以下にありますが、相当容量が大きいのでダウンロードするときは注意してください。
West Africa
https://gis1.servirglobal.net/TrainingMaterials/SAR/Data_ch3_wa.zip
HKH
https://gis1.servirglobal.net/TrainingMaterials/SAR/Data_ch3_hkh.zip

Switch to the data directory

os.chdir(datadirectory)
os.getcwd()  # Uncomment this line to display the present working directory
# glob.glob("*.vrt")   # Uncomment this line to see a List of the files 

os.chdirを使用して、作業ディレクトリを上記で設定したデータディレクトリに変更します。
os.getcwd()によって現在の作業ディレクトリを確認します。
コメントアウトされていますが、glob.glob("*.vrt")はvrtファイルのリスト表示します。

Acquisition Dates

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

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

dates=open(datefile).readlines()
tindex=pd.DatetimeIndex(dates)

# From the index we make and print a lookup table for 
# band numbers and 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()

日付ファイルを読み込み、Pandasの DatetimeIndex を使用して日付データをインデックス化します。
実施していることを簡単にリストすると以下を実施しています。

• 変数の初期化
  j を1で初期化しています。これはバンドのナンバリングに使用され、最初のバンドから順にナンバーが割り当てられます。
• ループの開始
  for i in tindex: は、tindex (DatetimeIndexオブジェクト) の各要素に対してループを行います。tindex は日付データのインデックスで、各バンドに対応する撮影日を含んでいます。
• 日付の表示
  print("{:4d} {}".format(j, i.date()), end=' ') は、バンド番号とその日付をフォーマットして表示します。{:4d} は整数を4桁で表示し、i.date() は日付データの日付部分のみを取り出します。end=' ' は、改行せずに次の出力を同じ行に続けるためのパラメータです。
• 改行の制御
  if j%5==1: print() は、5つのバンドごとに改行を挿入します。これにより、出力が整理されて見やすくなります。j を1ずつ増やすことでバンド番号を更新し、5の倍数+1になるたび（つまり、5バンド毎の後）改行が挿入されます。
  上記を実行すると以下のような結果になります。

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 

Image data

To open an image file and make it readable use the gdal.Open() function. This generates an image handle that can be used for further interactions with the file:

画像ファイルを開いて読めるようにするには、gdal.Open()関数を使います。この関数は特定の画像ファイルへの参照やポインタとして機能します。具体的には、画像ファイルを操作するために必要な情報や方法へのアクセスを提供するオブジェクトや変数を指します。一般的にそういった画像ファイルを操作するための情報をハンドルと呼びます。

img=gdal.Open(imagefile)

GDALライブラリを使用して、指定した画像ファイル（.vrt）を開きます。

To explore the image, e.g. number of bands, pixels, lines you can use several functions associated with the opened image object, e.g.:

画像のバンド数、ピクセル数、ライン数などを調べるには、以下のような関数が使用できます。

print(img.RasterCount) # Number of Bands
print(img.RasterXSize) # Number of Pixels
print(img.RasterYSize) # Number of Lines

画像のバンド数、ピクセル数、ライン数を表示します。結果は以下になります。

77
4243
3776

Reading data from an image band

To access any band in the image, use the img.GetRasterBand(x) function. E.g. to access the first band x=1, the last band: x=60.

画像内の任意のバンドにアクセスするには、img.GetRasterBand(x)関数を使用します。
例えば、最初のバンドにアクセスするにはx=1、最後のバンドにアクセスするにはx=60とします。

band=img.GetRasterBand(1)

第一バンドのデータを読み込みます。

Once a band is seleted, several functions associated with the band are available for further processing, e.g.
band.ReadAsArray(xoff=0,yoff=0,xsize=None,ysize=None)
So, to read the entire raster layer for the band:

一旦バンドが選択されると、そのバンドに関連するいくつかの機能が、さらなる処理のために利用可能となります。
例：band.ReadAsArray(xoff=0,yoff=0,xsize=None,ysize=None)
つまり、バンドのラスターレイヤー全体を読み込むには、次のようにします。

raster=band.ReadAsArray()

先ほど取得した第一バンドのデータを配列化しています。

Subsets

Because of the potentially large data volume when dealing with time series data stacks, it may be required to read only a subset of data.
With the gdal .ReadAsArray() function, subsets can be requested with offsets and size:
img.ReadAsArray(xoff=0, yoff=0, xsize=None, ysize=None)
xoff,yoff are the offsets from the upper left corner in pixel/line coordinates.
xsize,ysize specify the size of the subset in x-direction (left to right) and y-direction (top to bottom).
E.g., to read only a subset of 5x5 pixels with an offset of 5 pixels and 20 lines:

時系列データを扱う場合、データ量が膨大になる可能性があるため、データのサブセットのみを読み込みたい場合があります。
「gdal.ReadAsArray()」関数を以下のように使用すると、オフセットとサイズを指定してサブセットを作成できます。
img.ReadAsArray(xoff=0, yoff=0, xsize=None, ysize=None)
xoff,yoffは、左上隅からのピクセル/ライン座標でのオフセットです。
xsize,ysizeはX方向（左から右）とY方向（上から下）のサブセットのサイズを指定します。
例えば、5x5 ピクセルのサブセットだけを、5 ピクセルのオフセットと 20 行で読み込む場合は以下のように書きます。

オフセットとサイズは、画像データや配列データを部分的に取り扱う際に使用されるパラメータです。

オフセット（Offset）
オフセットは、特定の起点（通常は左上隅）からの相対的な距離を示します。画像やデータグリッドにおいてオフセットは、読み取りや処理を開始する特定の位置を指定するために使用されます。例えば、画像の左上隅から右に5ピクセル、下に20ピクセルの位置をオフセットとして指定することができます。この場合、xoff=5（横方向のオフセット）とyoff=20（縦方向のオフセット）となります。

サイズ（Size）
サイズは、データの読み取りや処理を行う領域の大きさを指定します。サイズは通常、幅（横方向の大きさ）と高さ（縦方向の大きさ）で定義されます。例えば、指定したオフセット位置から横方向に5ピクセル、縦方向に5ピクセルの範囲を読み取る場合、xsize=5、ysize=5と指定します。

例
オフセットとサイズを使用して特定のサブセットを読み取る例を考えます。例えば、画像データの中から、左上隅から右に5ピクセル、下に20ラインの位置を起点とし、そこから横5ピクセル、縦5ラインの範囲を抽出する場合、次のように指定します。
img.ReadAsArray(xoff=5, yoff=20, xsize=5, ysize=5)

raster_sub=band.ReadAsArray(5,20,5,5)

上記の説明の通り、band.ReadAsArray(xoff=5, yoff=20, xsize=5, ysize=5)を示しています。

The result is a two dimensional numpy array with the datatpye the data were stored in. We can inspect these data in python by simply typing the array name on the commandline:

結果は、2次元のnumpy配列になります。この配列はコマンドラインで配列名を入力するだけでデータ中身を表示することができます。

raster_sub

結果は以下になります。
この配列は5行5列の行列であり、画像の一部分のピクセル値を表しています。
各ピクセル値はバックスキャッター値を表しています。
バックスキャッターとは、合成開口レーダー（SAR）などのレーダー技術で使用される用語であり、レーダー波が物体や地表に当たって反射され、元の送信源（レーダー装置）に戻ってくる反射波のことを指します。

array([[4308, 4616, 4944, 4850, 4130],
       [3639, 4142, 4789, 5224, 4745],
       [3361, 3980, 4785, 5364, 4999],
       [3383, 3946, 4674, 5118, 4936],
       [3359, 3687, 4155, 4711, 5004]], dtype=uint16)

Displaying Bands in the Time Series of SAR Data

From the lookup table we know that bands 5 and 18 in the Niamey dataset are from late March and late October. Let's take look at these images.
HINT: Because python is an object oriented scripting language, we can often combine several steps (or function calls) into one command. See the trick below to access a raster band and read the data in one step.

上記のテーブル結果より、Niameyのデータセットのバンド5と18は3月下旬と10月下旬のものであることがわかります。これらの画像を見てみます。
HINT: pythonはオブジェクト指向のスクリプト言語なので、いくつかのステップ（または関数呼び出し）を1つのコマンドにまとめることができます。

lookup tableとは先述の「Acquisition Dates」項の結果のことを指しているのだと理解しましたが、結果を確認するに以下のようにバンド５は5月下旬であり、バンド18は10月下旬なので説明文が間違っているように見えます。

5 2015-05-21
18 2015-10-24
後述されるバンド5の結果表示でも2015-5-21と表示されているので多分誤記だと判断しスルーしました。
また、正直ここのヒントが何を指しているのかよくわかりませんでしたが、おそらく一般的なことしか言っていないのでここもスルーしました。

# These will select the two bands 
raster_1 = img.GetRasterBand(5).ReadAsArray()
raster_2 = img.GetRasterBand(18).ReadAsArray()

バンド5,18を配列データとして取得します。

Plotting in Python to Visualize the Image Bands

Matplotlib's plotting functions allow for powerful options to display imagery. We are following some standard approaches for setting up figures. First we are looking at a raster band and it's associated historgram.

Matplotlibのプロット機能は、画像を表示するための強力なオプションを提供します。ここでは、図の設定におけるいくつかの標準的なアプローチに従います。まず、ラスターバンドとそれに関連するヒストグラムを見ていきます。

fig = plt.figure(figsize=(16,8)) # Initialize figure with a size
ax1 = fig.add_subplot(121)  # 121 determines: 1 row, 2 plots, first plot
ax2 = fig.add_subplot(122)  # 122 determines: 1 row, 2 plots, second plot

# First plot: Image
bandnbr=5
ax1.imshow(raster_1,cmap='gray',vmin=2000,vmax=8000)
ax1.set_title('Image Band {} {}'.format(bandnbr,
                                    tindex[bandnbr-1].date()))

# Second plot: Historgram
# IMPORTANT: To get a histogram, we first need to *flatten* 
# the two-dimensional image into a one-dimensional vector.
h = ax2.hist(raster_1.flatten(),bins=100,range=(0,8000))
ax2.xaxis.set_label_text('Amplitude (Uncalibrated DN Values)')
_=ax2.set_title('Histogram Band {} {}'.format(bandnbr,
                                    tindex[bandnbr-1].date()))

個々のコードで実行していることを簡単にリストにするといかになります。

• 図の初期設定
  • fig = plt.figure(figsize=(16,8))
    16x8インチのサイズで新しい図（Figure）を初期化しています。
  • ax1 = fig.add_subplot(121) と ax2 = fig.add_subplot(122)
    1行2列のレイアウトで、それぞれ第1プロットと第2プロットを作成しています。
• 画像のプロット
  • bandnbr=5
    表示したい画像のバンド番号を指定しています。
  • ax1.imshow(raster_1, cmap='gray', vmin=2000, vmax=8000)
    raster_1 配列（バンド5の画像データ）をグレースケールで表示。
    vmin と vmax は画像の表示に用いるデータの最小値と最大値を設定しています。
  • ax1.set_title('Image Band {} {}'.format(bandnbr, tindex[bandnbr-1].date()))
    サブプロットのタイトルを設定しており、バンド番号とその日付を表示しています。

結果は以下のようになります。

ヒストグラムの形状が比較的中間、もしくはやや左寄りにあります。
植生はレーダー波を一定程度散乱させるため、中間の反射率を示すことが多いそうです。
若干左寄りにあるため、反射が弱い地域(平坦な地面や水面など)が若干多いかもしれません。

Writing a plotting function for the task

Below, the plotting commands used above are defined in a function named showImage. Several parameters can be passed to the function, some with default values listed at the end:
raster = a numpy two dimensional array
tindex = a panda index array for dates
bandnbr = the band number the corresponds to the raster
vmin = minimim value to display
vmax = maximum value to display
Note: By default, data will be linearly stretched between vmin and vmax.

上記で使用されたプロットようのコード(Plotting in Python to Visualize the Image Bands)を、「showImage」という名前の関数で定義します。
この関数には、以下の通り、いくつかのパラメータが渡され、いくつかはデフォルト値が設定されています。
raster: numpyの二次元配列
tindex: 日付のためのpandasインデックス配列
bandnbr: ラスターに対応するバンド番号
vmin: 表示する最小値
vmax: 表示する最大値
注記：デフォルトでは、データは vmin と vmax の間で線形に拡張されます。

def showImage(raster,tindex,bandnbr,vmin=None,vmax=None):
    fig = plt.figure(figsize=(16,8))
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)

    ax1.imshow(raster,cmap='gray',vmin=vmin,vmax=vmax)
    ax1.set_title('Image Band {} {}'.format(bandnbr,
                                tindex[bandnbr-1].date()))
    vmin=np.percentile(raster,2) if vmin==None else vmin #change vmin & vmax to change what values are displayed
    vmax=np.percentile(raster,98) if vmax==None else vmax
    ax1.xaxis.set_label_text(
        'Linear stretch Min={} Max={}'.format(vmin,vmax))
    
    
    h = ax2.hist(raster.flatten(),bins=100,range=(0,8000))
    ax2.xaxis.set_label_text('Amplitude (Uncalibrated DN Values)')
    ax2.set_title('Histogram Band {} {}'.format(bandnbr,
                                tindex[bandnbr-1].date()))

コードの解説は上記(Plotting in Python to Visualize the Image Bands)の通りなので割愛します。

EXERCISE 1: Read different bands and display them using the function showImage()

Use as a variable name for bands bandnbr. Use the already open image handle img to obtain the raster data from a band.

バンドの変数名としてbandnbrを使用してください。すでに開いている画像ハンドル img を使用して、バンドからラスターデータを取得します。
この指示は、プログラミングのコンテキストで、特定のバンド番号を指定して画像データを取り出す方法について説明しています。img は画像ファイルを開いた際に生成されたハンドルであり、このハンドルを通じて、bandnbr という変数で指定された特定のバンドのデータを抽出するプロセスを指しています。

# ENTER YOUR CODE HERE
showImage(raster_2,tindex,11,1000,2000)

上記で定義した「showImage」を用いてグラフを表示してみます。結果は以下のようになります。

EXERCISE 2: Read two different bands and display them side by side

The output should display two bands with a heading of the band numbers. Use the concept for figures with subplots from the function showImage(). Try your code to compare images from different years and different seasons.

バンド番号の見出し付きで二つのバンドを表示してみてください。showImage()関数からサブプロットを用いる概念を使用してください。異なる年や季節からの画像を比較するために、あなたのコードを試してみてください。

bandnbr1=51
raster1=img.GetRasterBand(bandnbr1).ReadAsArray()

bandnbr2=66
raster2=img.GetRasterBand(bandnbr2).ReadAsArray()

fig=plt.figure(figsize=(16,8))
ax1=fig.add_subplot(121)
ax2=fig.add_subplot(122)
ax1.imshow(raster1,vmin=2000,vmax=8000,cmap='gray')
ax2.imshow(raster2,vmin=2000,vmax=8000,cmap='gray')
ax1.set_title('Band {}   Date {}'.format(bandnbr1,tindex[bandnbr1-1].date()))
_=ax2.set_title('Band {}   Date {}'.format(bandnbr2,tindex[bandnbr2-1].date()))

Time Series Data Stacks

Just as we can use the ReadAsArray() function on a band, we can actually use it on the entire image data stack. To read an entire stack, i.e. all bands use the function on the image data handle:img.ReadAsArray()
CAUTION: Since this could potentially result in large memory need, it is wise to do some preliminary calcuations as to how large of a data set would be read in. for that we can do the following calculation:

For SAR data we typically use dataypes of
Float 32 bit (4 bytes per pixel) for power and dB data,
Unsigned Integer 16 bit (2 bytes per pixel) linearly scaled amplitudes, and
Unsigned Byte (1 byte per pixel) for dB-scaled to 8 bit data
The following table gives an overview of typically used data types for SAR data analysis in python:

Compare the result of the computation with the available RAM on the computer running the notebook.

バンドに対してReadAsArray()関数を使用するのと同様に、実際には画像データ全体に対してこの関数を使用することができます。データ全体、つまり全バンドを読み込むためには、画像データハンドルに対して「img.ReadAsArray()」関数を使用します。
注意：これにより大量のメモリが必要になる可能性があるため、読み込まれるデータセットのサイズを事前に計算することが好ましいです。そのためには、次の式にて容量を計算することができます。

データボリューム[GB] = バンド数 × ピクセル数 × ライン数 × バイト数 / 1024^3

SARデータでは通常、以下のデータタイプを使用します。

• 浮動小数点数32ビット（1ピクセルあたり4バイト）：パワーおよびdBデータ用
• 符号なし整数16ビット（1ピクセルあたり2バイト）：線形にスケールされた振幅用
• 符号なしバイト（1ピクセルあたり1バイト）：8ビットにスケールされたdBデータ用

以下の表は、PythonでのSARデータ解析に一般的に使用されるデータタイプの概要を示しています。

データタイプ	Numpy名	GDAL名	GDALコード	バイト数/ピクセル
浮動小数点数 32ビット	np.float32	gdal.GDT_Float32	6	4
符号なし整数 16ビット	np.uint16	gdal.GDT_UInt16	2	2
符号なしバイト	np.uint8	gdal.GDT_Byte	1	1

ノートブックを実行しているコンピュータの利用可能なRAMと計算結果を比較してみてください。

EXERCISE 3: Compute the Data Volume of the Raster Stack

Compute the estimated data volume from the data set opened with gdal.Open() using the img object information img.RasterXSize, img.RasterYSize, img.RasterCount, img.GetRasterBand(1).DataType

「gdal.Open()」で開いたデータセットから、imgオブジェクトの情報を使用して推定データボリュームを計算します。この計算には img.RasterXSize（画像の幅）、img.RasterYSize（画像の高さ）、img.RasterCount（画像に含まれるバンドの数）、img.GetRasterBand(1).DataType（最初のバンドのデータタイプ）を用います。
上記の式に従って、処理前に必要なメモリ量を確認します。

# Get the Data type
img.GetRasterBand(1).DataType

バンドのデータタイプを取得しています。符号なし16ビット整数が取得された場合のコードは以下になります。

#Use the lookup table for the number of bytes per pixel for this type:
bytespp=2
size=img.RasterCount*img.RasterXSize*img.RasterYSize*bytespp/(1024*1024*1024)
print('Data Volume for {}: {:.1f} Gigabytes'.format(img.GetDescription(),size))

Reading the SAR Time Series Subset

Let's read a image subset (offset 2000, 2000 / size 1024, 1024) of the entire time series data stack. The data are representatios of linearly scaled amplitudes scaled to unsigned 16 bit interger
We use the gdal ReadAsArray(xoff,yoff,xsize,ysize) function where
xoff = offset in pixels from upper left
yoff = offset in lines from upper left
xsize = number of pixels
ysize = number of lines
If ReadAsArray() is called without any parameters set, the entire image data stack is read. ReadAsArray() returns a numpy array of the form:

今まで説明されてきたことの繰り返しになります。
画像データの一部(オフセット 2000,2000 / サイズ 1024,1024)を読み込みます。
このデータは、符号なし16ビット整数にスケールされた振幅の線形表現です。
符号なし16ビット整数にスケールされた振幅の線形表現についてもう少し記述します。

• 符号なし16ビット整数 (Unsigned 16-bit Integer)
  これは、データが数値として表される形式の一つです。16ビット整数とは、各数値が16ビット（または2バイト）の情報を用いて表されることを意味します。符号なし（Unsigned）とは、負の数を含まず、0から65535までの正の数のみを表すことができる形式です。
• 振幅の線形表現
  振幅とは、一般的に波の強さや大きさを示す値です。例えば、今回の場合はバックスキャッターの振幅はその波の強さを表します。線形表現とは、振幅が実際の測定値と直線的（線形）な関係にあることを示します。これにより、データの解釈が直感的かつ簡単になります。
• スケールされたデータ
  データが「スケールされる」というのは、元の測定値をある範囲内で表現可能な形式に変換することを意味します。この場合、振幅が測定され、その値を0から65535の範囲にフィットするように調整（スケーリング）されます。このスケーリングによって、データはコンピュータやデジタルシステムで扱いやすくなります。

ここではgdalの「ReadAsArray(xoff, yoff, xsize, ysize)」関数を使います。パラメータについては説明済ですが、以下の通りです。
xoffは左上からのピクセル数でのオフセット(ずれ)です。
yoffは左上からのライン数でのオフセット(ずれ)です。
xsizeは読み込むピクセルの数です。
ysizeは読み込むラインの数です。
「ReadAsArray()」関数にパラメータを設定せずに呼び出すと、画像データの全スタックが読み込まれます。「ReadAsArray()」はnumpyの配列としてデータを返します。

# Alternatively you can make subset a tuple variable and use 
# it in the ReadAsArray function call prefixed with a star
subset=(2000,2000,1024,1024)
rasterDN = img.ReadAsArray(*subset)

まず「subset=(2000,2000,1024,1024)」にてsubsetという変数にタプルを割り当てています。このタプルは4つの要素を持ち、それぞれ xoff, yoff, xsize, ysize に対応します。
また「img.ReadAsArray()」はGDALライブラリの関数で、画像データをNumPy配列として読み込むために使われます。
(*subset)はPythonのアンパック演算子(*)を使っています。この演算子は、タプルやリストなどのシーケンス型の変数を関数の引数として展開する際に使用します。上記のコードでは、subset タプル内の各要素が「ReadAsArray()」関数の xoff, yoff, xsize, ysize という位置引数に直接渡されています。
変数subsetを調整するだけで異なる画像サブセットを簡単に読み込むことが可能になることを説明しています。

The numpy .shape object tells us the dimensions of this data stack as bands (here:time steps), lines, and pixels:

Numpyの「.shape」属性は、そのデータスタックの寸法をバンド（ここでは時間のステップ）、行、ピクセルとして教えてくれます。
時間のステップとは要は異なる連続した時点のことだと理解しています。

rasterDN.shape

rasterDNの寸法を確認すると、以下の結果が得られます。

(77, 1024, 1024)

時間のステップ、つまり時点としては77個の異なる時間のデータが存在することを示しています。
1024, 1024については、各画像のピクセルサイズを示しています。

Data conversion from linear scaled amplitudes to dB, power and amplitude data

The values of the raw image data show the linearly scaled amplitude values. These digital number (DN) values need to be converted to proper backscatter values of 𝛾𝑜 .
We consider conversion to dB scale (logarithmic scale) for the expression of the SAR backscatter, power, or amplitude scale.
SAR backscatter data of radiometrically terrain corrected data are often expressed as 𝜎𝑜 or the terrain flattened 𝛾𝑜 backscattering coefficients. For forest and land cover monitoring applications 𝛾𝑜 is the preferred metric.
Conversion from power to the logarithmic decibel (dB) scale follows:

As per widely used convention SAR backscatter data are often stored in 16bit unsigned integer values as linearly scaled amplitude data (referred to below as digital numbers DN), conversion to dB scale from the linear scaled amplitues is performed with a standard calibration factor of -83 dB. This is how ALOS SAR data are distributed by JAXA, how Earth Big Data LLC produces all SAR data including Sentinel-1 data, and how NISAR data will likely be scaled:
Conversion from amplitude to dB:

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

• デシベル（dB）スケールへの変換
  SARのバックスキャッターデータは、通常、デシベル（dB）という対数スケールで表現されます。これは、バックスキャッターの強さ（パワー）や振幅をより直感的に理解するのに役立ちます。
• バックスキャッター係数
  地形補正を施したSARデータは、𝜎𝑜や平坦化された 𝛾𝑜（バックスキャッター係数）で表現されることが多いです。森林や土地被覆のモニタリングでは、𝛾𝑜が指標として好まれます。
• dBスケールへの変換方法
  パワーからdBスケールへの変換は次の式で行います：

𝛾^𝑜_{dB} = 10 \times \log_{10}(𝛾𝑜_{power})

SARのバックスキャッターデータは、16ビット符号なし整数の線形スケール振幅データとして保存されることが多く、標準のキャリブレーションファクター「-83dB」を使って線形スケールからdBスケールへの変換が行われます。例えば、ALOS SARデータやEarth Big Data LLCが提供するSentinel-1データ、将来のNISARデータなどがこの方法で提供されます。

• 振幅からデシベルへの変換
  振幅からデシベルへの変換は以下の式で行います：

𝛾^𝑜_{dB} = 20 \times \log_{10}(DN) - 83

これらの変換を行うことで、データをより扱いやすくし、さまざまな環境条件下での変化を効果的に捉えることができます。

キャリブレーション係数「-83」について、この係数はセンサーによって変わるそうです。
そのセンサーの最低限界検出値が使われることが多いそうですが、なぜこの最低値をdBスケールでマイナスするのかがイマイチわかりませんでした。

デシベル値は次なような形式で計算されます。

\gamma^o_{dB} = 10 \times \log_{10}\left(\frac{P}{P_0}\right)

ここでPは測定したパワー、P_0は参照パワー（通常はセンサーの設計によって定義される）です。この式から、-83dBはセンサーが検出可能な最小の信号レベルを相対的に示しているそうですが、それが具体的に10^(-83)のパワーレベルを意味するわけではないようです。
あくまで-83dBが示すのは参照パワーP_0に対してどれだけ小さい信号がセンサーによって検出可能であるかの指標です。例えば、もしP_0を1ワットとした場合、-83dBの信号は、

P = 1 \times 10^{-8.3} [W]

と計算できるようです。
この参照パワーP_0は通常、センサーの設計仕様や校正プロセスに基づいて異なります。
説明文の計算式に参照パワーが表示されていないのは、既に絶対値として扱われているか、もしくは1という一般的な値が参照値として使用されているのなと思いましたが、真偽はわかりませんでした。
キャリブレーション係数が最低検出限界値の指標であるということはなんとなくわかりましたが、それをなぜマイナスするのかが分かりませんでした。この点について説明されたサイトや、回答をお持ちの方がいましたらご教示いただけると幸いです。
今はただセンサーによって決まったキャリブレーション係数を引いてスケールを行う、と単純記憶しておきます。

rasterdB=20*np.log10(rasterDN)-83

Conversion from dB to power:

dBスケールからパワースケールへの変換は以下の式で表せます。

\gamma^o_{\text{pwr}} = 10^{\frac{\gamma^o_{dB}}{10}}

rasterPwr=np.power(10.,rasterdB/10.)

Conversion from power to amplitude:

パワースケールから振幅への変換は以下の式で表せます。

\gamma^o_{\text{amp}} = \sqrt{\gamma^o_{\text{pwr}}}

rasterAmp=np.sqrt(rasterPwr)

Explore the image bands of the time steps

Let's explore how a band looks in the various image scales

バンドが様々なイメージスケールによってどのように見えるのか確認してみましょう。

Choose the band number and find which date it is

バンド番号を選び、それがどの日付であるかを調べます。

bandnbr=20
tindex[bandnbr-1]

tindexには日付データが入っています。バンド数から-1しているのは配列上、0から始まるためです。
結果は以下になります。

Timestamp('2015-11-29 00:00:00')

Below is the python code to create a four-part figure comparing the effect of the representation of the backscatter values in the DN, amplitude, power and dB scale.

以下は、DN、振幅、パワー、dBスケールでのバックスキャッター値の表現を比較する4つの図を作成するためのpythonコードです。

fig=plt.figure(figsize=(16,16))

ax1=fig.add_subplot(221)
ax2=fig.add_subplot(222)
ax3=fig.add_subplot(223)
ax4=fig.add_subplot(224)

ax1.imshow(rasterDN[bandnbr],cmap='gray',
           vmin=np.percentile(rasterDN,10),
           vmax=np.percentile(rasterDN,90))
ax2.imshow(rasterAmp[bandnbr],cmap='gray',
           vmin=np.percentile(rasterAmp,10),
           vmax=np.percentile(rasterAmp,90))
ax3.imshow(rasterPwr[bandnbr],cmap='gray',
           vmin=np.percentile(rasterPwr,10),
           vmax=np.percentile(rasterPwr,90))
ax4.imshow(rasterdB[bandnbr],cmap='gray',
           vmin=np.percentile(rasterdB,10),
           vmax=np.percentile(rasterdB,90))

ax1.set_title('DN Scaled (Amplitudes)')
ax2.set_title('Amplitude Scaled')
ax3.set_title('Power Scaled')
_=ax4.set_title('dB Scaled')

このコードは、Matplotlibを使用して、4つの異なる種類の画像データ（DNスケール、振幅スケール、パワースケール、デシベル（dB）スケール）を表示するためのものです。各画像はサブプロットに配置され、画像ごとに異なるスケールのデータを視覚的に比較できるようにしています。
図に4つのサブプロット（ax1, ax2, ax3, ax4）を追加します。これらは2x2のグリッド配置です（221, 222, 223, 224は位置を指定）。
各サブプロットは以下のグラフを表しています。

• ax1.imshow: DNスケールのデータ
• ax2.imshow: 振幅スケールのデータ
• ax3.imshow: パワースケールのデータ
• ax4.imshow: dBスケールのデータ

結果は以下のようになります。

Comparing histograms of the amplitude, power, and dB scaled data

それぞれのスケールを比較してみます。

# Setup for three part figure
fig=plt.figure(figsize=(16,4))
fig.suptitle('Comparison of Histograms of SAR Backscatter in Different Scales',fontsize=14)
ax1=fig.add_subplot(131)
ax2=fig.add_subplot(132)
ax3=fig.add_subplot(133)

# Important to "flatten" the 2D raster image to produce a historgram
ax1.hist(rasterAmp[bandnbr].flatten(),bins=100,range=(0.,0.6))
ax2.hist(rasterPwr[bandnbr].flatten(),bins=100,range=(0.,0.25))
ax3.hist(rasterdB[bandnbr].flatten(),bins=100,range=(-25,-5))

# Means, medians and stddev
amp_mean=rasterAmp[bandnbr].mean()
amp_std=rasterAmp[bandnbr].std()
pwr_mean=rasterPwr[bandnbr].mean()
pwr_std=rasterPwr[bandnbr].std()
dB_mean=rasterdB[bandnbr].mean()
dB_std=rasterdB[bandnbr].std()

# Some lines for mean and median
ax1.axvline(amp_mean,color='red')
ax1.axvline(np.median(rasterAmp[bandnbr]),color='blue')
ax2.axvline(pwr_mean,color='red',label='Mean')
ax2.axvline(np.median(rasterPwr[bandnbr]),color='blue',label='Median')
ax3.axvline(dB_mean,color='red')
ax3.axvline(np.median(rasterdB[bandnbr]),color='blue')

# Lines for 1 stddev
ax1.axvline(amp_mean-amp_std,color='gray')
ax1.axvline(amp_mean+amp_std,color='gray')
ax2.axvline(pwr_mean-pwr_std,color='gray',label='1 $\sigma$')
ax2.axvline(pwr_mean+pwr_std,color='gray')
ax3.axvline(dB_mean-dB_std,color='gray')
ax3.axvline(dB_mean+dB_std,color='gray')

ax1.set_title('Amplitude Scaled')
ax2.set_title('Power Scaled')
ax3.set_title('dB Scaled')
_=ax2.legend()

このコードは各スケールのヒストグラムを表示しています。処理内容をリストすると以下になります。

• fig=plt.figureで16x4インチの図を作成

• fig.suptitleで図全体のタイトルを設定

• fig.add_subplot(131)、fig.add_subplot(132)、fig.add_subplot(133)で3つのサブプロットを作成し、それぞれに異なるデータのヒストグラムをプロットするための領域を指定

• 各ヒストグラムは、対応するデータセット（振幅、パワー、dB）から選択されたバンドのデータを flattenメソッドで1次元配列に変換した後、histメソッドでプロット

  • flattenメソッドとは多次元の配列（例えば画像データは通常2次元です）を1次元配列に変換するためのメソッドです。この変換によって、画像データ内の各ピクセルが一列に並べられます。これを行う理由は、ヒストグラムを作成する際にデータが1次元のリストや配列である必要があるからです。
    たとえば、画像が1024x1024ピクセルのサイズであれば、その画像データは1024行と1024列で構成されています。flatten を使用すると、この1024x1024のデータが1,048,576個のピクセル値を持つ単一の長い列（1次元配列）に変換されます。
  • histメソッドは、数値データの配列を取り、そのデータの分布を表すヒストグラムを作成するための関数です。ヒストグラムはデータがどのように分布しているかを視覚的に示します。各棒（ビンと呼ばれます）は、特定の値の範囲にどれだけ多くのデータ点が存在するかを表します。

• bins=100はヒストグラムのビン（棒）の数を指定し、range でデータの表示範囲を設定

  • bins=100 とは、ヒストグラムを100個の棒で表示することを意味します。これにより、データの分布がより詳細に表示されます。
    range パラメータは、ヒストグラムに表示されるデータの最小値と最大値を設定します。例えば、range=(0., 0.6) は、ヒストグラムが0から0.6の範囲のデータ値のみを表示するようにします。

• 各データセットについて、平均、中央値、標準偏差を計算し、これらを縦線でヒストグラム上に表示

  • 赤線は平均、青線は中央値、灰色線は1標準偏差の範囲を表示

結果は以下のようになります。

Why is the scale important?

It is critical to use the correct scaling of SAR data for image processing operations. As we can see from the comparison of the histograms, the amplitude, power, and dB scales have different statistial distributions.
In time series analysis we often compare measurements at any given time step against the mean of the time series and compute its residuals. When we compute the mean of observations, it makes a difference whether we do that in power or dB scale. Since dB scale is a logarithmic scale, we cannot simply average data in that scale. Consider the following backscatter values and their mean:

Let's compute the mean of these values in power and dB scale and compare the result in dB scale:

スケールが重要である理由を説明します。SARデータの正しいスケーリングを使用することは、画像処理操作において非常に重要です。ヒストグラムの比較からわかるように、振幅、パワー、デシベル（dB）スケールはそれぞれ異なる統計的分布を示します。
時系列分析では、特定の時間での測定値を時系列の平均と比較し、その差異（残差）を計算することがよくあります。観測値の平均を計算する際には、それをパワースケールで行うか、dBスケールで行うかが重要です。dBスケールは対数スケールなので、そのスケールでデータを単純に平均化することはできません。以下の式は、バックスキャッター値とその平均の例です。
例を使って、これらの値の平均をパワーとdBスケールで計算し、結果をdBスケールで比較しましょう。
この説明は、異なるスケールでのデータ処理の違いがどのように結果に影響を与えるかを理解するためのものです。

g1_dB = -10
g2_dB = -15
g1_pwr = np.power(10.,-10/10.)
g2_pwr = np.power(10.,-15/10.)

mean_dB = (g1_dB+g2_dB)/2.
mean_pwr = (g1_pwr+g2_pwr)/2.
mean_pwr_inDB = 10.*np.log10(mean_pwr)

print('Mean averaging dB values         : {:.1f}'.format(mean_dB))
print('Mean averagin power values in dB : {:.1f}'.format(mean_pwr_inDB))

g1_dBとg2_dBはデシベルスケールでのバックスキャッター値です。
g1_pwrとg2_pwrは、デシベル値をパワースケールに変換したものです。
変換には

\gamma^o_1 = -10 \text{ dB}

\gamma^o_2 = -15 \text{ dB}

の公式を使用しています。
mean_dB はdBスケールの値を直接平均したものです。
mean_pwr はパワースケールでの値を平均したものです。
mean_pwr_inDB はパワースケールで平均化した値を再びデシベルスケールに変換しています。
結果は以下になります。

Mean averaging dB values         : -12.5
Mean averagin power values in dB : -11.8

dB値の平均(-12.5dB)は、二つのdB値そのままの平均です。
パワースケールでの平均化後のdB値(-11.8dB)は、それぞれのdB値をパワーに変換してから平均した後、再度dBに戻した値です。この値は直接平均したdB値よりも高いです。
この違いは、dBスケールが対数スケールであるため生じます。対数スケールでは、中間の値は直感的な算術平均よりも大きくなります。つまり、dB値を直接平均するのではなく、適切にパワースケールで平均を取り、それをdBに変換することが、より正確な平均を得る方法となります。
もう少し説明を加えると、dBスケールは対数的な変換を施したスケールであり、もとの振幅やパワー値がどのくらい大きいかを比較的感覚的に理解しやすくしてくれます。
しかし、その性質上、直接的な算術平均（各数値を足してその数で割る方法）を使用すると誤った結果を導くことがあります。具体的には、対数スケールでは値の幅が広がるため、低い数値と高い数値の間の平均が直感とは異なる結果を生み出します。
以下の引用も同じようなことを説明していますが、意訳も記述しておきます。

As one can see, there is a 0.7 dB difference in the average of these two 𝛾𝑜 backscatter values. If we make mean estimates of backscatter values, the correct scale in which operations need to be performed is the power scale. This is critical, e.g. when speckle filters are applied, spatial operations like block averaging are performed, or time series are analyzed. Very often we implement models that relate backscatter to biophysical variables like biomass, forest height, or use thresholds to determine change. Ensure that the proper scaling is done when working with the SAR data applying these models.
Another example of the effects can be illustrated with our backscatter data from the images we extracted. Consider a 1 hectare window extracted from our data sets with an off set of 500, 500 for band 20. We compute the mean over time and space of all the pixels.

この例では、2つのバックスキャッター値の平均値の違いが0.7dBあることがわかります。バックスキャッター値の平均推定を行う場合、操作を実行する正しいスケールはパワースケールです。これは、スペックルフィルターが適用されたり、ブロック平均のような空間操作が行われたり、時系列が分析される場合など、非常に重要です。また、バックスキャッターを生物物理的変数（例えば、バイオマス、森林の高さ）に関連付けるモデルを実装することがよくありますし、変化を判断するために閾値を使用することもあります。これらのモデルを適用する際には、SARデータの正しいスケーリングを確実に行うことが必要です。
難しい単語が出てきたので補足しておきます。

• スペックルフィルター（Speckle Filter）
  スペックルノイズは、レーダー画像やその他の画像タイプに見られる粒状のノイズパターンです。これは、画像の解像度によるものや、対象物の反射特性によるものなど、様々な要因で発生します。スペックルフィルターはこのノイズを軽減し、画像の質を向上させるために用いられます。
• ブロック平均（Block Averaging）
  ブロック平均は、画像やデータの特定のブロック（区画）にわたって平均値を計算し、データを単純化したり、ノイズを減らしたりする技術です。これにより、データの全体的な傾向がより明確になることがあります。
• 時系列分析（Time Series Analysis）
  時系列分析は、時間の経過とともに得られるデータ系列を分析する方法です。SARデータの場合、これを用いて地表や森林の変化を追跡したり、長期的な環境変動を監視したりします。
• バックスキャッターと生物物理的変数
  バックスキャッターは、レーダー波が地表や物体からどのように反射されるかを示す指標です。この値を生物物理的変数（例えば、バイオマスや森林の高さ）と関連付けるモデルは、環境の健康状態や変化を評価するために使用されます。
• 閾値（Threshold）
  変化検出やデータ分類のために設定される特定の値です。閾値を設定することで、バックスキャッター値がこの値を超えた場合にのみ、ある特定の処理を行うようにすることができます。

また、影響の例を説明するために、私たちが抽出した画像からのバックスキャッターデータを考えてみましょう。例として、データセットから500, 500のオフセットで1ヘクタールの画像を抽出し、その中の全ピクセルについて時空間の平均を計算します。
時空間の平均と言われてもよく分からないですが、時間と空間の両方の次元にわたってデータの平均値を計算することを指しているようです。具体的には、抽出した画像領域内の全ピクセル値について、時間を通じての平均（例えば、複数日または複数時刻にわたる画像データがある場合）と、空間的な平均（画像内の全ピクセルに対して）を同時に行っているようです。

offset=500
size=5
o1=offset
o2=offset+size

mean_dB = rasterdB[:,o1:o2,o1:o2].mean()
mean_dB

-11.302698

• rasterdBはdBスケールでのバックスキャッターデータを含む配列
• [ :, o1:o2, o1:o2]は全ての時間ステップにおいて、o1からo2までの範囲（正方形領域）のデータを抽出
• .mean()にて抽出したデータの平均値を計算

mean_pwr = rasterPwr[:,o1:o2,o1:o2].mean()
mean_pwr_in_dB = 10.* np.log10(mean_pwr)
mean_pwr_in_dB

-10.75519323348999

• rasterPwr はパワースケールでのデータを含む配列
• 上記同様に、特定の領域からデータを抽出し、その平均値を計算
• np.log10(mean_pwr)でパワー値をdBスケールに変換し、10倍

このコードは、特定の範囲内の画像データからバックスキャッターの平均dBスケール値と平均パワースケール値を計算し、それをdBに変換する処理を行っています。
この結果から、直接dB値で計算した平均と、パワー値をdBに変換してから計算した平均にはわずかながら違いが見られます。
この違いは、前述の説明の通り、dBスケールの対数的性質によるものです。パワースケールで計算した後にdBスケールへ変換することにより、より正確な平均値を取得できます。

As one can see, a difference of more than 0.5 dB is found simply by operating in the different scales. Hence: CAUTION!
何度も説明されていますが、異なるスケールで計算を行うだけで、0.5dB以上の差が生じることがわかります。
ですので、注意が必要とのことです。

Exploring Polarization Differences

We look at the backscatter characteristics in SAR data from like-polarized (same transmit and receive polarzation, hh or vv) and cross-polarized (vh or hv polarization). For this, we read a timestep in both polarizations, plot the histograms, and display the images in dB scale. First, we open the images, pick the bands from the same acquisition date, read the raster bands and convert them to dB scale.

SARデータにおける偏光の違いを調査します。ここでは、同じ偏光（送信と受信が同じ偏光、例えばhhやvv）と異なる偏光（vhやhv偏光）のバックスキャッターの特性を見ていきます。このために、両方の偏光にて、特定時点のデータを読み込み、ヒストグラムをプロットし、dBスケールで画像を表示します。まず、画像を開いて、同じ取得日のバンドを選び、ラスターバンドを読み込んでdBスケールに変換します。

# Open the Images
img_like=gdal.Open(imagefile)
img_cross=gdal.Open(imagefile_cross)
# Pick the bands, read rasters and convert to dB
bandnbr_like=20
bandnbr_cross=20
rl=img_like.GetRasterBand(bandnbr_like).ReadAsArray()
rc=img_cross.GetRasterBand(bandnbr_cross).ReadAsArray()
rl_dB=20.*np.log10(rl)-83
rc_dB=20.*np.log10(rc)-83

このコードは、SAR画像データの同一偏光と異なる偏光のバンドを開き、それらを読み込んでdBスケールに変換する手順を示しています。それぞれのステップを以下にリストします。

• gdal.Open関数を使用して、指定されたファイルパスの画像ファイルを取得
• img_likeは同一偏光（例えばhhまたはvv）の画像を保持
• img_crossは異なる偏光（例えばvhまたはhv）の画像を保持
• 各画像から特定のバンド（この例では20番バンド）を選択
• GetRasterBandメソッドは、指定されたバンド番号に対応するバンドオブジェクトを取得
• ReadAsArray は、そのバンドのデータをNumPy配列に変換
• rlは同一偏光の画像データ、rcは異なる偏光の画像データを保持
• データをdBスケールに変換(dB変換には一般的に 20 * log10(データ) の式を使用)
• -83は標準のキャリブレーションファクターで、実際の値は使用するデータやセンサーによって異なる
• rl_dBとrc_dBはそれぞれ同一偏光と異なる偏光のデータをdBスケールに変換した結果を保持

Now, we explore the differences in the polarizations by plotting the images with their histograms. We look at the dB ranges over which the histograms spread, and can adjust the linear scaling in the image display accordingly to enhace contrast. In the case below
C-vv like-polarized data are mostly spread from -17.5 to -5 dB
C-vh cross-polarized data are mostly spread from -25 to -10 dB
Thus, we note that the cross-polarized data exhibit a larger dynamic range of about 2.5 dB

異なる偏光を持つSARデータの違いを探るために、それぞれの画像とそれに対応するヒストグラムをプロットして比較します。ヒストグラムを用いて、どのdB範囲にデータが集中しているかを見ます。その情報を基に、画像のコントラストを強化するための調整を行います。
具体的なデータ範囲としては、同じ偏光のデータ（C-vv）は主に-17.5dBから-5dBの範囲に、異なる偏光のデータ（C-vh）は-25dBから-10dBの範囲に広がっています。異なる偏光のデータはより広い範囲（約2.5dBの違い）を持っており、これによりデータから得られる情報の幅が広がっています。

fig,ax=plt.subplots(nrows=2,ncols=2,figsize=(16,16))
fig.suptitle('Comaprison of Like- and Cross-Polarized Sentinel-1 C-band Data',
             fontsize=14)
ax[0][0].set_title('C-VV Image')
ax[0][1].set_title('C-VH Image')
ax[1][0].set_title('C-VV Histogram')
ax[1][1].set_title('C-VH Histogram')
ax[0][0].axis('off')
ax[0][1].axis('off')
ax[0][0].imshow(rl_dB,vmin=-17.5,vmax=-5,cmap='gray')
ax[0][1].imshow(rc_dB,vmin=-25,vmax=-10,cmap='gray')
ax[1][0].hist(rl_dB.flatten(),range=(-25,-5),bins=100)
ax[1][1].hist(rc_dB.flatten(),range=(-25,-5),bins=100)
fig.tight_layout()  # Use the tight layout to make the figure more compact

このPythonコードは、Matplotlibライブラリを使用して、Sentinel-1のC-bandデータの同一偏光（C-VV）と異なる偏光（C-VH）の画像、およびそれらのヒストグラムを表示するためのものです。以下に各部分の説明を行います。

• plt.subplots
  2x2のグリッドにサブプロットを作成しています。このグリッドは、上の行に画像を、下の行にヒストグラムを配置するためのものです。
• figsize=(16,16)
  フィギュアのサイズを16インチ四方に設定しています。
• fig.suptitle
  フィギュア全体のタイトルを設定しています。
• ax[0][0].set_title('C-VV Image') ~ ax[1][1].set_title('C-VH Histogram')
  各サブプロットにタイトルを設定しています。
• ax[0][0].axis('off')
  axis('off') で不要な軸の表示をオフにしています。
• imshow
  画像を表示しており、vminとvmaxを設定することで表示の際のカラースケールの範囲を定義しています。これにより画像のコントラストを強調できます。
• flattenで画像データを1次元配列に変換後、histでヒストグラムを描画
• tight_layout
  サブプロット間のスペースを自動調整し、見た目を整えています。

結果は以下になります。

左上の画像（C-VV Image）と右上の画像（C-VH Image）を比較すると、テクスチャやコントラストに若干の違いが見られます。C-VVの画像は少し明るく見えるのに対し、C-VHの画像はやや暗く、コントラストが高いようです。
これらの違いは、C-VVは通常、表面の水平方向の特性をより詳細に捉え、C-VHは地表の垂直方向の特性や異質性をより強調します。
左下のヒストグラム（C-VV Histogram）と右下のヒストグラム（C-VH Histogram）は、それぞれの画像のピクセル値の分布を示しています。
両方のヒストグラムは鋭いピークを持っており、データが特定のdB範囲に集中していることを示しています。これは、データが特定の反射特性を持つ地表の特徴に対応していることを意味します。
C-VVのヒストグラムはやや右にシフトしており、全体的に高い反射率を示しています。これに対して、C-VHのヒストグラムは左に広がりがあり、より広いダイナミックレンジを示しています。これは、C-VHが異なるタイプの地表特性やより複雑な地形条件を捉えていることを示唆しています。
つまるところ、以下のようなことが言えるようです。

• C-VV偏光（同一偏光）
  一般的に地表の滑らかな部分や均一な特性を反映します。画像が比較的明るく、反射が強い地域をはっきりと見せるため、水面や平坦な土地など、滑らかな特性を持つ地表の検出に適しています。
• C-VH偏光（異なる偏光）
  地表の粗さや複雑な構造をより詳細に捉えることができます。これは、画像が全体的に暗く、ヒストグラムの広がりが大きいことからもわかります。つまり、植生の密度が高い地域や都市部など、様々な方向からの反射が複雑に絡み合う地域の特性をより詳しく捉えることができるます。
• 画像解析の適用性
  これらの違いを理解することで、特定の偏光データを然るべき用途で使用することが可能になります。例えば、水資源の管理や農地の監視にはC-VVデータが、森林のバイオマス推定や都市開発のモニタリングにはC-VHデータが適していると考えられます。
• データのダイナミックレンジ
  C-VHデータが示す広いダイナミックレンジは、そのデータが持つ情報の豊富さを示しています。広い範囲のdB値は、異なる種類の地表が混在する地域を詳しく解析する際に有用です。

EXERCISE 4: Explore different Seasons in different polarizations

Change the band numbers bandnbr_like and bandnbr_cross in the cell above to explore different bands.

まとめ

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