警告: このチュートリアルで説明する Landsat ETM+ データと OLI データを調和させる手順は古く、Landsat Collection 2 地表反射率データを扱う場合は推奨されず、必要ありません。詳細については、Landsat の調和に関するよくある質問をご覧ください。コレクション 1 のデータは非推奨となり、このチュートリアルのコードは機能しなくなり、更新もされなくなります。
このチュートリアルでは、Landsat ETM+ の地表反射率を Landsat OLI の地表反射率に調和させる方法について説明します。内容は以下のとおりです。
- スペクトル変換関数
- 分析に対応したデータを作成する関数
- 時系列の可視化の例
このガイドは、35 年以上にわたる Landsat データの地域時系列を調和させて可視化するためのエンドツーエンドのガイドです。このガイドで説明する手法は、対象の地域にすぐに適用できます。
ETM+ から OLI への変換係数は TM にも適用されます。そのため、このチュートリアルでは、ETM+ への参照は TM と同義です。
Landsat について
Landsat は、1972 年から中程度の解像度の地球画像を収集している衛星画像プログラムです。最長の宇宙ベースの地球観測プログラムとして、景観変化の時空間的傾向を特定するための貴重な時間記録を提供します。このチュートリアルでは、Thematic Mapper（TM）、Enhanced Thematic Mapper Plus（ETM+）、Operational Land Imager（OLI）の機器データを使用します。これらは密接に関連しており、比較的簡単に統合して、1984 年から現在までの連続した記録を生成する一貫した時系列にすることができます。この時系列は、センサーごとに 16 日の頻度で、30 メートルの空間解像度で生成されます。マルチスペクトル スキャナ機器は、Landsat の記録を 1972 年まで遡りますが、このチュートリアルでは使用しません。データが大きく異なるため、後続のセンサーとの統合が困難です。すべてのセンサー間のハーモナイゼーションの例については、Savage 他（2018）と Vogeler 他（2018）を参照してください。
ハーモナイゼーションの理由
Roy ら（2016 年）は、用途によっては、Landsat ETM+ と OLI のスペクトル特性にわずかながらも有意な差があることを示しています。データセットを調和させる理由としては、Landsat TM、ETM+、OLI にまたがる長い時系列の作成、ETM+ SLC-off ギャップや雲/影のマスキングによる欠測値の影響を軽減するためのほぼ同日の年内コンポジットの生成、時系列内の観測頻度の増加などがあります。詳しくは、上記のリンク先の原稿をご覧ください。
手順
関数
以下は、ETM+ を OLI に調和させ、このチュートリアルのアプリケーション セクションで使用してピクセルのスペクトル時系列履歴を可視化する分析対応データを生成するために必要な一連の関数です。
ハーモナイゼーション
調和は、Roy et al.（2016）の表 2 の OLS 回帰係数に従って、ETM+ スペクトル空間から OLI スペクトル空間への線形変換によって実現されます。バンドごとの係数は、次のディクショナリで、傾き（slopes）と切片（itcps）の画像定数として定義されています。なお、y 切片の値は USGS Landsat 地表反射率データのスケーリングに合わせて 10,000 倍されます。
var coefficients = {
itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172])
.multiply(10000),
slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
同じスペクトル応答ウィンドウの ETM+ と OLI のバンド名が等しくないため、標準化する必要があります。次の関数は、各データセットから反射率バンドと pixel_qa バンドのみを選択し、それらが表す波長範囲に従って名前を変更します。
// Function to get and rename bands of interest from OLI.
function renameOli(img) {
return img.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
// Function to get and rename bands of interest from ETM+.
function renameEtm(img) {
return img.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
最後に、変換関数を定義します。この関数は、線形モデルを ETM+ データに適用し、OLI との整合性を保つためにデータ型を Int16 としてキャストし、後で雲と影のマスキングで使用するために pixel_qa バンドを再アタッチします。
function etmToOli(img) {
return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
.multiply(coefficients.slopes)
.add(coefficients.itcps)
.round()
.toShort()
.addBands(img.select('pixel_qa'));
}
雲と影のマスキング
分析準備完了データでは、雲と雲の影がマスクされている必要があります。次の関数は、CFmask（Zhu et al.、2015）
各 Landsat USGS 地表反射率画像に含まれる pixel_qa バンド。雲と雲の影として識別されたピクセルを null に設定します。
function fmask(img) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = img.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask)
.eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return img.updateMask(mask);
}
スペクトル指標の計算
このアプリケーションでは、正規化バーン比率（NBR）スペクトル指標を使用して、山火事の影響を受けた森林のピクセルのスペクトル履歴を表します。NBR が使用されるのは、Cohen et al.（2018）が、13 個のスペクトル指標/バンドの中で、米国全土の森林破壊（シグナル）に関して NBR の信号対雑音比が最も高いことを発見したためです。次の関数で計算されます。
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
独自のアプリケーションでは、別のスペクトル インデックスを使用することもできます。追加のインデックスを変更または追加する場所は次のとおりです。
関数を組み合わせる
各データセットのラッパー関数を定義して、上記のすべての関数を統合し、それぞれの画像コレクションに適用しやすくします。
// Define function to prepare OLI images.
function prepOli(img) {
var orig = img;
img = renameOli(img);
img = fmask(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
// Define function to prepare ETM+ images.
function prepEtm(img) {
var orig = img;
img = renameEtm(img);
img = fmask(img);
img = etmToOli(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
アプリケーションでは、関数を含めるか除外するかを決定できます。必要に応じてこれらの関数を変更します。
期間グラフの例
上記の prepOli 関数と prepEtm 関数は、Landsat 地表反射率コレクションにマッピングして、クロスセンサー分析対応データを作成し、ピクセルまたはピクセル領域のスペクトル時系列を可視化できます。この例では、35 年以上の時系列を作成し、単一ピクセルのスペクトル履歴を表示します。この特定のピクセルは、1980 年代と 1990 年代に何らかの摂動を受け、2011 年に大規模な火災が発生した、成熟した太平洋岸北西部の針葉樹林のパッチ（図 1）の最近の履歴に関連しています。
図 1. 関心領域の例の場所とサイトの文字。米国オレゴン州フッド山の北斜面にある、成熟した太平洋岸北西部の針葉樹林。画像提供: Google Earth、米国農務省森林局、Landsat、Copernicus。
対象領域（AOI）を定義する
このアプリケーションの結果は、ピクセルの Landsat 観測の時系列です。ee.Geometry.Point オブジェクトは、ピクセルの位置を定義するために使用されます。
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
アプリでは、別のピクセルを選択する場合があります。上記の経度と緯度の座標を置き換えてください。または、ee.Geometry.Polygon() や ee.Geometry.Rectangle() などの他の ee.Geometry オブジェクト定義を使用して、ピクセル グループのスペクトル履歴を要約することもできます。詳しくは、デベロッパー ガイドのジオメトリのセクションをご覧ください。
Landsat センサー コレクションを取得する
OLI、ETM+、TM の Landsat USGS 地表反射率コレクションを取得します。
各データセットの詳細については、リンクをご覧ください。
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
イメージ コレクション フィルタを定義する
フィルタを定義して、対象領域の空間境界、光合成のピーク シーズン、品質で画像コレクションを制限します。
var colFilter = ee.Filter.and(
ee.Filter.bounds(aoi), ee.Filter.calendarRange(182, 244, 'day_of_year'),
ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
ee.Filter.or(
ee.Filter.eq('IMAGE_QUALITY', 9),
ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));
これらのフィルタ条件は必要に応じて変更してください。上記のデータの説明ページの [画像プロパティ] タブで、フィルタする他のプロパティを確認します。
コレクションを準備する
コレクションを統合し、フィルタを適用して、すべての画像に prepImg 関数をマッピングします。次のスニペットの結果は、フィルタ条件を満たし、分析可能な NBR に処理された OLI、ETM+、TM センサー コレクションの画像を含む単一の ee.ImageCollection です。
// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOli);
etmCol = etmCol.filter(colFilter).map(prepEtm);
tmCol = tmCol.filter(colFilter).map(prepEtm);
// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);
すべての観測値を表示する期間グラフを作成する
センサー間の調和は、すべての観測値を散布図としてプロットし、センサーを色で区別することで、定性的に評価できます。Earth Engine の ui.Chart.feature.groups 関数は、このユーティリティを提供しますが、まず観測された NBR 値をプロパティとして各画像に追加する必要があります。AOI と交差するすべてのピクセルのメジアンを計算し、結果を画像プロパティとして設定する画像コレクションにリージョン縮小関数をマッピングします。
var allObs = col.map(function(img) {
var obs = img.reduceRegion(
{geometry: aoi, reducer: ee.Reducer.median(), scale: 30});
return img.set('NBR', obs.get('NBR'));
});
独自のアプリケーションで、AOI が広い場合は、scale を増やし、reduceRegion 関数の bestEffort、maxPixels、tileScale 引数を指定して、オペレーションが最大ピクセル数、メモリ、タイムアウトの制限を超えないようにすることを検討してください。中央値の削減関数を、任意の統計情報に置き換えることもできます。詳細については、デベロッパー ガイドの画像領域の統計情報のセクションをご覧ください。
コレクションを ui.Chart.feature.groups 関数で受け取れるようになりました。この関数には、コレクションと、グラフにマッピングするコレクション オブジェクトのプロパティの名前が必要です。ここで、'system:time_start' プロパティは x 軸変数として、'NBR' は y 軸変数として機能します。NBR プロパティは、縮小される画像のバンド名に基づいて、上記の reduceRegion 関数によって命名されています。別の帯域（「NBR」以外）を使用する場合は、それに応じて y 軸のプロパティ名引数を変更します。最後に、グループ化（系列）変数を「SATELLITE」に設定します。これには個別の色が割り当てられます。
var chartAllObs =
ui.Chart.feature.groups(allObs, 'system:time_start', 'NBR', 'SATELLITE')
.setChartType('ScatterChart')
.setSeriesNames(['TM', 'ETM+', 'OLI'])
.setOptions({
title: 'All Observations',
colors: ['f8766d', '00ba38', '619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
pointSize: 6,
dataOpacity: 0.5
});
print(chartAllObs);
処理が完了すると、図 2 のようなグラフがコンソールに表示されます。注意事項:
- 1 年に複数回の観測が行われます。
- 時系列は 3 つの異なるセンサーで構成されています。
- センサーバイアスは検出されませんでした。
- 観測頻度は年によって異なります。
- 年内変動はありますが、年ごとに一貫しており、比較的狭い範囲に収まっています。NDVI で試してみると（img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');）、年内および年間の応答の変動がはるかに大きくなります。
- NBR は、ほとんどの時系列で高い値を示し、わずかな摂動が見られます。
- 森林火災が原因で NBR の反応が大幅に減少していることが、
- ただし、火災（ダラー湖）は 2011 年 9 月に発生しました。検出年のエラーは、年間の複合期間が 7 月から 8 月であるためです。この期間後に発生した変更は、次に利用可能な null 以外の観測値（翌年以降になる可能性があります）まで検出されません。
- 植生の回復（NBR レスポンスの増加）は、NBR の大幅な損失が検出されてから 2 年後に始まります。
図 2. Landsat TM、ETM+、OLI の表現を含む単一ピクセルのスペクトル応答時系列グラフ。TM 画像と ETM+ 画像は、線形変換によって OLI に調和されます。画像は 7 月から 8 月のもので、高画質でフィルタされています。
年間の中央値を表示する期間グラフを作成する
時系列を簡素化してノイズを除去するには、年内の観測値の削減を適用できます。ここでは中央値が使用されています。
まず、画像を年ごとにグループ化します。各画像の ee.Image.Date からコレクション設定の「year」をマッピングして、各画像に新しい「year」プロパティを追加します。
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
新しい「year」プロパティを使用すると、同じ年の画像を結合でグループ化できます。個別の画像年コレクションと完全なコレクションの結合により、同じ年のリストが提供され、そこから中央値の削減を実行できます。完全なコレクションを、個別の年の代表のセットにサブセット化します。
var distinctYearCol = col.distinct('year');
個別の年コレクション（distinctYearCol）と完全なコレクション（col）の間で一致する年を特定するフィルタと結合を定義します。一致は「year_matches」というプロパティに保存されます。
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
結合を適用し、結果の FeatureCollection を ImageCollection に変換します。
var joinCol = ee.ImageCollection(join.apply(distinctYearCol, col, filter));
// Apply median reduction among matching year collections.
var medianComp = joinCol.map(function(img) {
var yearCol = ee.ImageCollection.fromImages(img.get('year_matches'));
return yearCol.reduce(ee.Reducer.median())
.set('system:time_start', img.date().update({month: 8, day: 1}));
});
特定の年の画像のコレクションを表す統計情報として中央値以外のものを使用する場合は、上記の ee.Reducer を変更します。
最後に、年内夏の観測セットから NBR の中央値を表示するグラフを作成します。必要に応じてリージョン削減統計情報を変更します。
var chartMedianComp = ui.Chart.image
.series({
imageCollection: medianComp,
region: aoi,
reducer: ee.Reducer.median(),
scale: 30,
xProperty: 'system:time_start',
})
.setSeriesNames(['NBR Median'])
.setOptions({
title: 'Intra-annual Median',
colors: ['619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
lineWidth: 6
});
print(chartMedianComp);
処理が完了すると、図 3 のようなグラフがコンソールに表示されます。
図 3. Landsat TM、ETM+、OLI のデータで表された単一ピクセルの夏のスペクトル応答時間の中央値の時系列グラフ。TM 画像と ETM+ 画像は、線形変換によって OLI に調和されます。画像は 7 月から 8 月のもので、高品質なものに絞り込まれています。
その他の情報
代替変換関数
Roy et al.（2016）の表 2 には、ETM+ の地表反射率を OLI の地表反射率に変換する、またはその逆を行うための OLS 回帰係数と RMA 回帰係数が示されています。上記のチュートリアルでは、OLS 回帰による ETM+ から OLI への変換のみを示しています。すべての翻訳オプションの関数は、コード エディタ スクリプトまたは GitHub ソース スクリプトにあります。
ダラー レイクの火災
このチュートリアルで言及している火災は、2011 年 9 月中旬に米国オレゴン州フッド山の北斜面で発生したダラー湖火災です。詳しくは、こちらのブログ投稿と Varner et al.（2015）をご覧ください。