公告：凡是在 2025 年 4 月 15 日前註冊使用 Earth Engine 的非商業專案，都必須驗證非商業用途資格，才能繼續存取。如未在 2025 年 9 月 26 日前完成驗證，存取權可能會暫停。

本頁面由 Cloud Translation API 翻譯而成。

Landsat ETM+ 到 OLI 的調和

作者： jdbcode

警告：本教學課程中介紹的 Landsat ETM+ 和 OLI 資料協調程序已過時，不建議使用，且使用 Landsat Collection 2 地表反射率資料時也不需要。詳情請參閱有關 Landsat 協調的常見問題。「集合 1」資料已淘汰，本教學課程中的程式碼將無法運作或更新。

在程式碼編輯器中開啟

本教學課程的重點是將 Landsat ETM+ 地表反射率調和為 Landsat OLI 地表反射率。以下是這個平台的特點：

頻譜轉換函式
可建立分析就緒資料的函式
時間序列視覺化範例

這份指南旨在提供端對端指引，協助您協調及視覺化 35 年以上的 Landsat 區域時間序列資料，並立即套用至感興趣的區域。

請注意，將 ETM+ 轉換為 OLI 的係數也適用於 TM。因此，在本教學課程中，提及 ETM+ 時，即等同於 TM。

關於 Landsat

Landsat 是一項衛星成像計畫，自 1972 年以來持續收集地球中等解析度圖像。Landsat 是歷史最悠久的太空地球觀測計畫，可提供寶貴的時間記錄，用於識別地景變化的時空趨勢。本教學課程會使用主題繪圖機 (TM)、增強型主題繪圖機 Plus (ETM+) 和作業陸地成像儀 (OLI) 的儀器資料。這兩項資料集密切相關，且相對容易整合成一致的時間序列，可產生 1984 年至今的連續記錄，每個感應器每 16 天提供一次資料，空間解析度為 30 公尺。多光譜掃描儀可將 Landsat 記錄追溯至 1972 年，但本教學課程不會使用這項儀器。其資料與後續感應器差異甚大，因此整合起來相當困難。如需所有感應器之間相互協調的範例，請參閱 Savage 等人 (2018 年) 和 Vogeler 等人 (2018 年) 的研究。

為什麼要進行協調

Roy 等人 (2016 年) 證明，視應用程式而定，Landsat ETM+ 和 OLI 的光譜特徵之間存在微小但可能顯著的差異。您可能需要協調資料集的原因包括：產生涵蓋 Landsat TM、ETM+ 和 OLI 的長期時序、產生近期的年內複合資料，以減少 ETM+ SLC 關閉間隙和雲/陰影遮罩造成的觀測資料遺失影響，或提高時序內的觀測頻率。如需更多資訊，請參閱上方的連結手稿。

操作說明

函式

以下是一系列函式，可將 ETM+ 調整為 OLI，並產生可用於分析的資料，以便在本教學課程的應用程式部分中，以視覺化方式呈現象素的光譜時間歷程記錄。

統合

根據 Roy 等人 (2016 年) 表 2 OLS 迴歸係數，將 ETM+ 光譜空間線性轉換為 OLI 光譜空間，即可完成調和作業。以下列字典定義各頻帶的係數，並使用斜率 (slopes) 和截距 (itcps) 圖像常數。請注意，y 截距值會乘以 10,000，以符合 USGS Landsat 地表反射率資料的縮放比例。

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+ 資料、將資料類型轉換為 Int16，以確保與 OLI 一致，並重新附加 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 等人，2015)，並隨附於每張 Landsat USGS 地表反射率影像，可將識別為雲朵和雲影的像素設為空值。pixel_qa

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 等人 (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 年以上的時間序列，並顯示單一像素的光譜記錄。這個特定像素與太平洋西北部成熟針葉樹林地塊的近期歷史有關 (圖 1)，該地塊在 1980 年代和 1990 年代經歷了一些擾動，並在 2011 年發生大規模燃燒。

感興趣的區域

圖 1. 位置和網站字元，例如感興趣的區域。美國奧勒岡州胡德山北坡的成熟太平洋西北針葉林。圖片來源：Google 地球、美國農業部林務局、Landsat 和 Copernicus。

定義感興趣的區域 (AOI)

這項應用程式的結果是像素的 Landsat 觀測時間序列。ee.Geometry.Point 物件用於定義像素的位置。

var aoi = ee.Geometry.Point([-121.70938, 45.43185]);

您可能會在應用程式中選取其他像素。方法是替換上述經緯度座標。或者，您可以使用其他 ee.Geometry 物件定義 (例如 ee.Geometry.Polygon() 和 ee.Geometry.Rectangle())，匯總一組像素的光譜記錄。詳情請參閱開發人員指南的「幾何圖形」一節。

擷取 Landsat 感應器集合

取得 OLI、ETM+ 和 TM 的 Landsat USGS 地表反射率集合。OLI TM

點選連結即可進一步瞭解各個資料集。

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 函式對應至所有圖片。下列程式碼片段的結果是單一 ee.ImageCollection，其中包含符合篩選條件的 OLI、ETM+ 和 TM 感應器集合中的圖片，並經過處理，可供分析 NBR。

// 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 軸變數。請注意，上述 reduceRegion 函式是根據要縮減的圖片中的頻帶名稱，命名 NBR 屬性。如果使用其他頻帶 (而非「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 的圖表。請注意以下幾點：

每年會進行多次觀測。
時間序列由三種不同的感應器組成。
沒有可辨識的感應器偏差。
每年觀察的頻率不一。
每年都有一些年內差異，但差異相對較小。使用 NDVI 嘗試此操作 (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) 會顯示更大的年內和年際反應變異性。
在大多數時間序列中，NBR 仍維持在高點，僅出現些微擾動。
森林大火造成 NBR 回應大幅減少，
但請注意，火災 (多拉湖) 發生於 2011 年 9 月，偵測年份發生錯誤，是因為年度綜合日期範圍為 7 月至 8 月；如果變更發生在這個時間範圍之後，系統要等到下一個可用的非空值觀測結果 (可能是下一年或更晚) 才會偵測到變更。
偵測到 NBR 大幅減少後，植被會在兩年後開始復原 (NBR 回應增加)。

時間序列所有觀測項目

圖 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);