衛星データ（衛星画像）のDN（デジタル・ナンバー）を反射率に変換する方法

反射率（大気上端）

反射率（大気上端）への変換方法

　衛星データのDN（デジタル・ナンバー）は、画像として見るには良いですが、画像を解析するにはDNを反射率に変換する必要があります。
　DNを大気上端反射率に変換するには、DNを輝度に変換する式と、輝度を反射率に変換する式を使います。


１．輝度への変換

【ALOSの場合、Terra/ASTER、GeoEye1の場合】
　以下にDNを輝度へ変換する式を示します。


　L = a・DN＋b


　ここで、
　L：輝度（W/m2/sr/μm）・・・ALOS、Terra/ASTER
　L：輝度（W/cm2/sr/μm）・・・GeoEye1
　DN：デジタル・ナンバー
　a：ゲイン
　b：オフセット


　式から分かるように、DNを輝度に変換するにはゲインとオフセットを知る必要があります。
　ゲインとオフセットは、ALOS衛星の場合は、衛星データを入手したときのリーダファイルのアンシラリ２に記載されているため、ここから情報を得ることになります。
　ASTERの場合は、拡張子が「.vnirsh」のテキストファイルに記載されています。


【WorldView-2の場合】

　DN（qのこと）を輝度に変換する式は、次のようになります。


　L = absCalFactor・q／Δλ


　ここで、
　L：輝度（W/m2/sr/μm）
　q：デジタル・ナンバー
　absCalFactor：補正係数
　Δλ：有効バンド幅


　absCalFactorと、Δλは、データを入手したときにテキストファイルに記載されています。


【WorldView-3の場合】

　DNを輝度に変換する式は、次のようになります。


　L = GAIN・DN・（absCalFactor／Δλ）＋OFFSET


　ここで、
　L：輝度（W/m2/sr/μm）
　DN：デジタル・ナンバー
　absCalFactor：補正係数
　Δλ：有効バンド幅
　GAIN：乗数
　OFFSET：加数


　absCalFactorと、Δλは、データを入手したときにテキストファイルに記載されています。


　GAINは、以下のとおりです。


　Pnchromatic：0.923
　Coastal：0.863
　Blue：0.905
　Green：0.907
　Yellow：0.938
　Red：0.945
　Red Edge：0.980
　NIR1：0.982
　NIR2：0.954
　SWIR1：1.160
　SWIR2：1.184
　SWIR3：1.173
　SWIR4：1.187
　SWIR5：1.286
　SWIR6：1.336
　SWIR7：1.340
　SWIR8：1.392


　OFFSETは、以下のとおりです。


　Pnchromatic：-1.700
　Coastal：-7.154
　Blue：-4.189
　Green：-3.287
　Yellow：-1.816
　Red：-1.350
　Red Edge：-2.617
　NIR1：-3.752
　NIR2：-1.507
　SWIR1：-4.479
　SWIR2：-2.248
　SWIR3：-1.806
　SWIR4：-1.507
　SWIR5：-0.622
　SWIR6：-0.605
　SWIR7：-0.423
　SWIR8：-0.302


２．大気上端反射率への変換

　輝度を反射率に変換する式は、次のようになります。


　ρ = (L・π・d2）／（Esun・cosθ）


　ここで、
　ρ：反射率（%）
　L：輝度（W/m2/sr/μm）
　d：地球太陽間距離（AU）
　Esun：分光太陽照度（W/m2/μm）
　θ：太陽天頂角



　また、Esunは、衛星およびセンサによって異なります。



LandsatのEsunは、USGSのサイト(https://landsat.usgs.gov/esun)より抜粋しています。

Landsat 1-5 MSS
　バンド1: 1848
　バンド2: 1588
　バンド3: 1235
　バンド4: 856.6

LANDSAT 4 TM
　バンド1: 1958
　バンド2: 1826
　バンド3: 1554
　バンド4: 1033
　バンド5: 214.7
　バンド7: 80.70

LANDSAT 5 TM
　バンド1: 1958
　バンド2: 1827
　バンド3: 1551
　バンド4: 1036
　バンド5: 214.9
　バンド7: 80.65

LANDSAT 7 ETM+
　バンド1: 1970
　バンド2: 1842
　バンド3: 1547
　バンド4: 1044
　バンド5: 225.7
　バンド7: 82.06
　バンド8: 1369

ALOS AVNIR-2
　バンド１：1943.3
　バンド２：1813.7
　バンド３：1562.3
　バンド４：1076.5

ASTER
　バンド1：1848
　バンド2：1549
　バンド3：1114
　バンド4：225.4

GeoEye1
　Panchromatic:161.7 (mW/cm2/μm) 　（単位に注意）
　Blue：196.0
　Green:185.3
　Red：150.5
　Near IR:103.9

　上記の数値を他の衛星と同じ単位にすると、
　Panchromatic:1617 (W/m2/μm)
　Blue：1960
　Green:1853
　Red：1505
　Near IR:1039

WorldView-2
　Panchromatic：1580.8140 (W/m2/μm)
　Coastal:1758.2229
　Blue：1974.2416
　Green:1856.4104
　Yellow:1738.4791
　Red：1559.4555
　Red Edge：1342.0695
　NIR1：1069.7302
　NIR2：861.2866

WorldView-3
【WRCの場合】
　Pnchromatic：1583.58 (W/m2/μm)
　Coastal：1743.81
　Blue：1971.48
　Green：1856.26
　Yellow：1749.4
　Red：1555.11
　Red Edge：1343.95
　NIR1：1071.98
　NIR2：863.296
　SWIR1：494.595
　SWIR2：261.494
　SWIR3：230.518
　SWIR4：196.766
　SWIR5：80.365
　SWIR6：74.7211
　SWIR7：69.043
　SWIR8：59.8224


　この式から分かるように、輝度から反射率を求めるには、地球太陽間距離dと太陽天頂角θを知る必要があります。

　太陽天頂角は、ALOSの場合は、リーダファイルのシーンヘッダに記載されている太陽仰角（太陽高度）、あるいはサマリ情報ファイル（summary.txt）に記載されている太陽仰角から、以下の式により求めることができます。


　太陽天頂角 = 90度　−　太陽仰角


　また、地球太陽間距離d[単位：AU]は、イメージ撮影日(世界標準時）から以下の式で、Julian Dayを求めます。


　year＝観測年（西暦）、month＝観測月、day＝観測日、UT＝観測時＝hour＋minute/60＋second/3600

　ただし、monthが、1月または2月の場合、year＝観測年-1、month＝観測月＋12 を用いる。

　Julian Day（JD）は、上記の　year、month、day、UTと、下記のAおよびBより、次のように求められます。

　A＝int(year/100)　

　B＝2−A＋int(A/4)

　JD = int{365.25☓（year＋4716）｝＋int｛30.6001☓（month＋１）｝＋day＋UT/24＋B−1524.5

　この値が、紀元前4713年1月1日の正午からの日数（JD）になります。

　「US. Naval Observatory」によると、このJDから地球太陽間距離は、次式で求めることができます。

　D＝JD−2451545.0

　g＝357.529＋0.98560028☓D　・・・・単位：度

　gをラジアン単位にするには、g＝g☓π／180度

　よって、地球太陽間距離d [AU] は、

　d＝1.00014−0.01671☓cos（g）−0.00014☓cos（2☓g）


　地球太陽間距離のエクセルの計算ファイルは、ここからダウンロードできます。


　当サイトでは、ALOS衛星AVNIR-2のリーダファイルからAVNIR-2反射率計算パラメータ表示するフリーソフトを提供しています。

　なお、当サイトで紹介しているリモートセンシング画像処理ソフトRSPは、32ビット（float）RAWファイル画像を取り扱う（上記の演算処理ができるとともに作成した画像を表示する）ことができます。RSPにおいて、ビットマップファイルをRAW（float）ファイルに変換する方法は、メニュー「RSP」→「RSP Q＆A」の「放射輝度、反射率」を参照してください。


Landsat-8における反射率への変換方法

　Landsat-8では、16ビットのQcal（デジタル・ナンバー）を以下の式を使って、反射率に変換します。


　ρλ = （Ｍρ・Ｑcal＋Ａρ）／cos(90°−θSE）


　ここで、ρλ：反射率（%）

　　　　　Ｑcal：デジタル・ナンバー

　　　　　Ｍρ：乗数

　　　　　Ａρ：加数

　　　　　θSE：太陽高度(度）


　詳しくは、Using the USGS Landsat 8 Productを見てください。


　各値は、データをダウンロードしたときのテキストファイルに記載されています。


　Landsat-8を例として、RSP を使った反射率変換方法を以下に記載します。


　具体の反射率への変換方法は、Landsat-8画像を見てみように記載しているため、参考にしてください。


-----------------------------------------------------------


【変換式を使って求める場合】

　1.TIFFファイルをRAW（16ビット）ファイルに変換する

　　メニュー「ファイル」→「フォーマット変換」→「TIFF(GeoTIFF）」→「16bit→RAW（unsigned int）」でRAWファイルに変換します。


　2.RAW(unsigned int)ファイルをRAW(float）ファイルに変換する

　　メニュー「ファイル」→「フォーマット変換」→「RAW」→「RAW（unsigned int)→RAW(float）」で、変換します。


　3.Qcalを反射率に変換する

　　メニュー「ツール」→「float処理」→「演算」→「＝a*(A)+ b」で、aの値にＭρを、bの値にＡρを入力します。

　　次に、メニュー「ツール」→「float処理」→「演算」→「＝a*(A)+ b」で、aの値に（1÷sin（θSE））を、bの値に0（ゼロ）を入力します。


　4.ゼロ以下、１以上の値を修正する

　　反射率の範囲は、理論上、0から1までの数値です。しかし、上記の変換式で、ゼロ以下の値や、1以上の値が求まる場合があります。

　　このため、ゼロ以下の値をゼロに、1以上の値を1に補正しなければなりません。

　　この補正処理を「ピクセル値変更」処理を使って行います。

　　メニュー「ツール」→「float処理」→「編集」→「ピクセル値変更」で、マイナス領域（範囲）のデータを0にします。

　　続けて、メニュー「ツール」→「float処理」→「編集」→「ピクセル値変更」で、1以上のデータを1にします。


　以上の処理を経て、Landsat-8のデジタル・ナンバーを反射率に変換することができます。



-----------------------------------------------------------


【RSPの専用機能で求める場合】


　1.TIFFファイルをRAW（16ビット）ファイルに変換する

　　メニュー「ファイル」→「フォーマット変換」→「TIFF(GeoTIFF）」→「16bit→RAW（Unsigned int）」でRAWファイルに変換します。


　2.RAW(unsigned int)ファイルを反射率画像（RAW，Float）ファイルに変換する

　　メニュー「ファイル」→「ツール」→「Landsat-8（RAW，unsigned int)→反射率」で、Ｍρ、Ａρ、太陽高度（太陽仰角）（度）の値を入力し、出力先として「Float」を選択し、「作成」ボタンを押します。

　　「1.」で変換したRAWファイルを選択し、次に、保存先ファイル名（反射率画像ファイル名）を入力します。


　これにより、Landsat-8のデジタル・ナンバーを反射率に変換することができます。なお、「2.」の処理は繰り返し行うことができます。

　また、反射率画像（RAW，Unsigned int）ファイルに変換する場合は、上記「2.」の出力先選択として「Unsigned integer (*10000)」を選択します。

　「*10000」は、反射率（最大値：1）を10000倍した値で出力する（最大値：10000）という意味です。