Pengumuman: Semua project nonkomersial yang terdaftar untuk menggunakan Earth Engine sebelum 15 April 2025 harus memverifikasi kelayakan nonkomersial untuk mempertahankan akses. Jika Anda belum melakukan verifikasi hingga 26 September 2025, akses Anda mungkin ditangguhkan.

Halaman ini diterjemahkan oleh Cloud Translation API.

Penyelarasan Landsat ETM+ ke OLI

Edit di GitHub

Laporkan masalah

Histori halaman

Penulis: jdbcode

PERINGATAN: Prosedur yang dijelaskan dalam tutorial ini untuk menyelaraskan data Landsat ETM+ dan OLI sudah tidak berlaku dan tidak direkomendasikan atau diperlukan saat bekerja dengan data reflektansi permukaan Landsat Collection 2. Untuk mengetahui informasi selengkapnya, lihat entri FAQ tentang harmonisasi Landsat. Data Collection 1 tidak digunakan lagi dan kode dalam tutorial ini tidak akan berfungsi atau diupdate lagi.

Buka Di Editor Kode

Tutorial ini membahas penyelarasan pantulan permukaan Landsat ETM+ dengan pantulan permukaan Landsat OLI. Referensi ini menyertakan:

Fungsi transformasi spektral
Fungsi untuk membuat data yang siap dianalisis
Contoh visualisasi deret waktu

Panduan ini dimaksudkan sebagai panduan end-to-end untuk menyelaraskan dan memvisualisasikan deret waktu regional data Landsat selama lebih dari 35 tahun yang dapat langsung diterapkan ke wilayah yang Anda minati.

Perhatikan bahwa koefisien untuk mengubah ETM+ menjadi OLI juga berlaku untuk TM. Oleh karena itu, dalam tutorial ini, referensi ke ETM+ identik dengan TM.

Tentang Landsat

Landsat adalah program pencitraan satelit yang telah mengumpulkan citra Bumi beresolusi sedang sejak tahun 1972. Sebagai program pengamatan Bumi berbasis ruang angkasa terlama, program ini memberikan catatan temporal yang berharga untuk mengidentifikasi tren spatiotemporal dalam perubahan lanskap. Data instrumen Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), dan Operational Land Imager (OLI) digunakan dalam tutorial ini. Data ini terkait erat dan relatif mudah digabungkan menjadi deret waktu yang konsisten yang menghasilkan catatan berkelanjutan dari tahun 1984 hingga saat ini, dengan irama 16 hari per sensor dengan resolusi spasial 30 meter. Instrumen Multispectral Scanner memperpanjang catatan Landsat hingga tahun 1972, tetapi tidak digunakan dalam tutorial ini. Datanya cukup berbeda, sehingga integrasi dengan sensor berikutnya menjadi sulit. Lihat Savage et al. (2018) dan Vogeler et al. (2018) untuk mengetahui contoh harmonisasi di semua sensor.

Alasan harmonisasi

Roy et al. (2016) menunjukkan bahwa ada perbedaan kecil, tetapi berpotensi signifikan antara karakteristik spektral Landsat ETM+ dan OLI, bergantung pada penerapannya. Alasan Anda mungkin ingin menyelaraskan set data meliputi: membuat deret waktu panjang yang mencakup Landsat TM, ETM+, dan OLI, membuat komposit intra-tahunan yang mendekati tanggal untuk mengurangi efek pengamatan yang hilang dari kesenjangan SLC-off ETM+ dan masking awan/bayangan, atau meningkatkan frekuensi pengamatan dalam deret waktu. Lihat naskah yang ditautkan di atas untuk mengetahui informasi selengkapnya.

Petunjuk

Fungsi

Berikut adalah serangkaian fungsi yang diperlukan untuk menyelaraskan ETM+ dengan OLI dan menghasilkan data siap analisis yang akan digunakan di bagian penerapan tutorial ini untuk memvisualisasikan histori spektral-temporal piksel.

Harmonisasi

Harmonisasi dicapai melalui transformasi linear ruang spektral ETM+ ke ruang spektral OLI sesuai dengan koefisien yang disajikan dalam Roy et al. (2016) Tabel 2 koefisien regresi OLS. Koefisien khusus band ditentukan dalam kamus berikut dengan konstanta gambar kemiringan (slopes) dan intersep (itcps). Perhatikan bahwa nilai intersep y dikalikan dengan 10.000 agar sesuai dengan penskalaan data reflektansi permukaan Landsat USGS.

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])
};

Nama band ETM+ dan OLI untuk jendela respons spektral yang sama tidak sama dan perlu distandardisasi. Fungsi berikut hanya memilih band reflektansi dan band pixel_qa dari setiap set data, lalu mengganti namanya sesuai dengan rentang panjang gelombang yang diwakilinya.

// 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']);
}

Terakhir, tentukan fungsi transformasi, yang menerapkan model linear ke data ETM+, mentransmisikan jenis data sebagai Int16 agar konsisten dengan OLI, dan melampirkan kembali band pixel_qa untuk digunakan nanti dalam masking awan dan bayangan.

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'));
}

Masking awan dan bayangan

Data yang siap dianalisis harus memiliki awan dan bayangan awan yang ditutupi. Fungsi berikut menggunakan CFmask (Zhu et al., 2015) pixel_qa disertakan dengan setiap gambar reflektansi permukaan USGS Landsat untuk menyetel piksel yang diidentifikasi sebagai awan dan bayangan awan ke 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);
}

Penghitungan indeks spektral

Aplikasi yang akan datang menggunakan indeks spektral rasio pembakaran yang dinormalisasi (NBR) untuk merepresentasikan histori spektral piksel hutan yang terkena dampak kebakaran hutan. NBR digunakan karena Cohen et al. (2018) menemukan bahwa di antara 13 indeks/band spektral, NBR memiliki rasio sinyal terhadap derau yang paling besar terkait gangguan hutan (sinyal) di seluruh AS. Hal ini dihitung dengan fungsi berikut.

function calcNbr(img) {
  return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}

Dalam aplikasi Anda sendiri, Anda dapat memilih untuk menggunakan indeks spektral yang berbeda. Di sini adalah tempat untuk mengubah atau menambahkan indeks tambahan.

Menggabungkan fungsi

Tentukan fungsi wrapper untuk setiap set data yang menggabungkan semua fungsi di atas agar mudah diterapkan ke koleksi gambar masing-masing.

// 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()));
}

Dalam aplikasi, Anda dapat memutuskan untuk menyertakan atau mengecualikan fungsi. Ubah fungsi ini sesuai kebutuhan.

Contoh rangkaian waktu

Fungsi wrapper prepOli dan prepEtm di atas dapat dipetakan ke koleksi pantulan permukaan Landsat untuk membuat data siap analisis lintas sensor guna memvisualisasikan kronologi spektral piksel atau wilayah piksel. Dalam contoh ini, Anda akan membuat deret waktu lebih dari 35 tahun dan menampilkan histori spektral untuk satu piksel. Piksel khusus ini terkait dengan sejarah terbaru petak hutan konifer dewasa di Pacific Northwest (Gambar 1) yang mengalami gangguan pada tahun 1980-an dan 1990-an serta kebakaran dengan intensitas tinggi pada tahun 2011.

Area minat

Gambar 1. Karakter lokasi dan situs untuk contoh area yang diminati. Hutan konifer dewasa di Pacific Northwest di lereng utara Mt Hood, OR, Amerika Serikat. Gambar atas izin: Google Earth, USDA Forest Service, Landsat, dan Copernicus.

Menentukan area minat (AOI)

Hasil aplikasi ini adalah deret waktu pengamatan Landsat untuk piksel. Objek ee.Geometry.Point digunakan untuk menentukan lokasi piksel.

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

Dalam aplikasi, Anda mungkin memilih piksel yang berbeda. Lakukan dengan mengganti koordinat bujur dan lintang di atas. Atau, Anda dapat meringkas histori spektral sekelompok piksel menggunakan definisi objek ee.Geometry lainnya, seperti ee.Geometry.Polygon() dan ee.Geometry.Rectangle(). Lihat bagian Geometri di Panduan Developer untuk mengetahui informasi selengkapnya.

Mengambil koleksi sensor Landsat

Dapatkan koleksi pantulan permukaan USGS Landsat untuk OLI, ETM+, dan TM.

Buka link untuk mempelajari setiap set data lebih lanjut.

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

Menentukan filter koleksi gambar

Tentukan filter untuk membatasi kumpulan gambar berdasarkan batas spasial area yang diminati, musim fotosintesis puncak, dan kualitas.

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

Ubah kriteria pemfilteran ini sesuai keinginan Anda. Identifikasi properti lain yang akan difilter berdasarkan tab "Properti Gambar" di halaman deskripsi data yang ditautkan di atas.

Menyiapkan koleksi

Gabungkan koleksi, terapkan filter, dan petakan fungsi prepImg ke semua gambar. Hasil cuplikan berikut adalah satu ee.ImageCollection yang berisi gambar dari koleksi sensor OLI, ETM+, dan TM yang memenuhi kriteria filter, dan diproses ke NBR yang siap dianalisis.

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

Buat diagram deret waktu yang menampilkan semua pengamatan

Kesesuaian antar-sensor dapat dinilai secara kualitatif dengan memetakan semua pengamatan sebagai diagram sebar yang membedakan sensor berdasarkan warna. Fungsi ui.Chart.feature.groups Earth Engine menyediakan utilitas ini, tetapi pertama-tama nilai NBR yang diamati harus ditambahkan ke setiap gambar sebagai properti. Memetakan fungsi reduce region di atas kumpulan gambar yang menghitung median semua piksel yang berpotongan dengan AOI dan menetapkan hasilnya sebagai properti gambar.

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'));
});

Di aplikasi Anda sendiri, jika AOI Anda adalah area yang luas, Anda dapat mempertimbangkan untuk meningkatkan scale dan menentukan argumen bestEffort, maxPixels, tileScale untuk fungsi reduceRegion guna memastikan bahwa operasi tidak melebihi batasan piksel, memori, atau waktu tunggu maksimum. Anda juga dapat mengganti pengurang median dengan statistik pilihan Anda. Lihat bagian Statistik Wilayah Gambar di Panduan Developer untuk mengetahui informasi selengkapnya.

Koleksi kini dapat diterima oleh fungsi ui.Chart.feature.groups. Fungsi ini mengharapkan koleksi dan nama properti objek koleksi untuk dipetakan ke diagram. Di sini, properti 'system:time_start' berfungsi sebagai variabel sumbu x dan 'NBR' sebagai variabel sumbu y. Perhatikan bahwa properti NBR diberi nama oleh fungsi reduceRegion di atas berdasarkan nama band dalam gambar yang dikurangi. Jika Anda menggunakan band lain (bukan 'NBR'), ubah argumen nama properti sumbu y dengan tepat. Terakhir, tetapkan variabel pengelompokan (seri) ke 'SATELLITE', yang memiliki warna berbeda.

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

Diagram yang mirip dengan Gambar 2 akan muncul di konsol setelah beberapa waktu pemrosesan. Perhatikan beberapa hal:

Ada beberapa pengamatan per tahun.
Deret waktu terdiri dari tiga sensor yang berbeda.
Tidak ada bias sensor yang dapat dibedakan.
Frekuensi pengamatan bervariasi setiap tahun.
Ada beberapa variasi dalam setahun, yang konsisten di seluruh tahun dan relatif kecil. Mencoba ini dengan NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) menunjukkan variabilitas respons intra- dan antar-tahunan yang jauh lebih besar.
NBR tetap tinggi untuk sebagian besar deret waktu dengan beberapa gangguan kecil.
Penurunan respons NBR yang signifikan dan cepat akibat kebakaran hutan terlihat pada
Namun, perhatikan bahwa kebakaran (Dollar Lake) terjadi pada September 2011. Error pada tahun deteksi disebabkan oleh rentang tanggal gabungan tahunan dari bulan Juli hingga Agustus; perubahan yang terjadi setelah periode ini tidak akan diambil hingga pengamatan non-null berikutnya tersedia, yang bisa terjadi pada tahun berikutnya atau lebih lama.
Pemulihan vegetasi (peningkatan respons NBR) dimulai dua tahun setelah kehilangan NBR besar terdeteksi.

Deret waktu semua observasi

Gambar 2. Diagram deret waktu respons spektral untuk satu piksel dengan representasi dari Landsat TM, ETM+, dan OLI. Gambar TM dan ETM+ diselaraskan dengan OLI melalui transformasi linier. Gambar diambil dari bulan Juli hingga Agustus dan difilter untuk kualitas tinggi.

Buat diagram deret waktu yang menampilkan median tahunan

Untuk menyederhanakan dan menghilangkan gangguan dari deret waktu, pengurangan pengamatan dalam tahunan dapat diterapkan. Di sini, median digunakan.

Langkah pertama adalah mengelompokkan gambar menurut tahun. Tambahkan properti 'year' baru ke setiap gambar dengan memetakan 'year' setelan koleksi dari ee.Image.Date setiap gambar.

var col = col.map(function(img) {
  return img.set('year', img.date().get('year'));
});

Properti 'year' baru dapat digunakan untuk mengelompokkan gambar dari tahun yang sama dengan gabungan. Penggabungan antara kumpulan tahun gambar yang berbeda dan kumpulan lengkap memberikan daftar tahun yang sama, yang dapat digunakan untuk melakukan pengurangan median. Buat subkumpulan koleksi lengkap menjadi serangkaian representasi tahun yang berbeda.

var distinctYearCol = col.distinct('year');

Tentukan filter dan gabungan yang akan mengidentifikasi tahun yang cocok antara kumpulan tahun unik (distinctYearCol) dan kumpulan lengkap (col). Kecocokan akan disimpan ke properti bernama 'year_matches'.

var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');

Terapkan gabungan dan konversi FeatureCollection yang dihasilkan ke 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}));
});

Ubah ee.Reducer di atas jika Anda lebih memilih statistik selain median untuk mewakili kumpulan gambar dalam tahun tertentu.

Terakhir, buat diagram yang menampilkan nilai NBR median dari kumpulan pengamatan musim panas dalam setahun. Ubah statistik pengurangan region sesuai keinginan Anda.

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

Diagram yang mirip dengan Gambar 3 akan muncul di konsol setelah beberapa waktu pemrosesan.

Median deret waktu

Gambar 3. Diagram deret waktu respons spektral musim panas median untuk satu piksel dengan representasi dari Landsat TM, ETM+, dan OLI. Gambar TM dan ETM+ diselaraskan dengan OLI melalui transformasi linier. Gambar diambil dari bulan Juli hingga Agustus dan difilter untuk mendapatkan kualitas tinggi.

Informasi tambahan

Fungsi transformasi alternatif

Tabel 2 Roy et al. (2016) memberikan koefisien regresi OLS dan RMA untuk mengubah pantulan permukaan ETM+ menjadi pantulan permukaan OLI dan sebaliknya. Tutorial di atas hanya menunjukkan transformasi ETM+ ke OLI dengan regresi OLS. Fungsi untuk semua opsi terjemahan dapat ditemukan di skrip Editor Kode atau skrip sumber GitHub.

Kebakaran Dollar Lake

Kebakaran yang disebutkan dalam tutorial ini adalah kebakaran Dollar Lake, yang terjadi pada pertengahan September 2011 di lereng utara Mt Hood, OR, AS. Buka postingan blog ini dan lihat Varner et al. (2015) untuk mengetahui informasi selengkapnya tentang hal ini.