Regresi dengan Kumpulan Data Sematan Satelit

Edit di GitHub
Laporkan masalah
Histori halaman
Tutorial ini adalah bagian dari serangkaian tutorial tentang set data Sematan Satelit, lihat juga Pengantar, Klasifikasi Tanpa Pengawasan, Klasifikasi dengan Pengawasan, dan Penelusuran Kemiripan.

Kolom penyematan dapat digunakan sebagai input/prediktor fitur untuk regresi dengan cara yang sama seperti yang digunakan untuk klasifikasi.

Dalam tutorial ini, kita akan mempelajari cara menggunakan lapisan kolom sematan 64D sebagai input untuk analisis regresi berganda yang memprediksi biomassa di atas permukaan tanah (AGB).

Misi Global Ecosystem Dynamics Investigation (GEDI) NASA mengumpulkan pengukuran LIDAR di sepanjang transek darat pada resolusi spasial 30 m dengan interval 60 m. Kita akan menggunakan set data GEDI L4A Raster Aboveground Biomass Density yang berisi estimasi titik kepadatan biomassa di atas permukaan tanah (AGBD) yang akan digunakan sebagai variabel yang diprediksi dalam model regresi.

Pilih region

Mari kita mulai dengan menentukan region yang diminati. Untuk tutorial ini, kita akan memilih region di Western Ghats, India, dan menentukan poligon sebagai variabel geometri. Atau, Anda dapat menggunakan Alat Gambar di Editor Kode untuk menggambar poligon di sekitar wilayah yang diinginkan yang akan disimpan sebagai variabel geometri di Impor. Kami juga menggunakan peta dasar Satelit, yang memudahkan Anda menemukan area yang ditumbuhi vegetasi.

var geometry = ee.Geometry.Polygon([[
  [74.322, 14.981],
  [74.322, 14.765],
  [74.648, 14.765],
  [74.648, 14.980]
]]);

// Use the satellite basemap
Map.setOptions('SATELLITE');


Gambar: Memilih area yang diminati untuk prediksi biomassa di atas permukaan tanah

Pilih jangka waktu

Pilih tahun yang ingin kita jalankan regresinya. Perlu diingat bahwa Sematan Satelit diagregasi pada interval tahunan, jadi kita menentukan periode untuk sepanjang tahun.

var startDate = ee.Date.fromYMD(2022, 1, 1);
var endDate = startDate.advance(1, 'year');

Menyiapkan set data Embedding Satelit

Gambar Sematan Satelit 64-band akan digunakan sebagai prediktor untuk regresi. Kita memuat set data Sematan Satelit, memfilter gambar untuk tahun dan wilayah yang dipilih.

var embeddings = ee.ImageCollection('GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL');

var embeddingsFiltered = embeddings
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

Gambar Sematan Satelit dikelompokkan ke dalam petak dan ditayangkan dalam proyeksi untuk zona UTM petak. Hasilnya, kita mendapatkan beberapa petak Sematan Satelit yang mencakup wilayah yang diminati. Untuk mendapatkan satu gambar, kita perlu menggabungkannya. Di Earth Engine, gabungan gambar input diberi proyeksi default, yaitu WGS84 dengan skala 1 derajat. Karena kita akan menggabungkan dan memproyeksikan ulang mozaik ini nanti dalam tutorial, sebaiknya pertahankan proyeksi aslinya. Kita dapat mengekstrak informasi proyeksi dari salah satu petak dan menyetelnya pada mosaik menggunakan fungsi setDefaultProjection().

// Extract the projection of the first band of the first image
var embeddingsProjection = ee.Image(embeddingsFiltered.first()).select(0).projection();

// Set the projection of the mosaic to the extracted projection
var embeddingsImage = embeddingsFiltered.mosaic()
  .setDefaultProjection(embeddingsProjection);

Menyiapkan mozaik GEDI L4A

Karena perkiraan biomassa GEDI akan digunakan untuk melatih model regresi kami, sangat penting untuk memfilter data GEDI yang tidak valid atau tidak andal sebelum menggunakannya. Kami menerapkan beberapa masker untuk menghapus pengukuran yang berpotensi salah.

  • Hapus semua pengukuran yang tidak memenuhi persyaratan kualitas (l4_quality_flag = 0 dan degrade_flag > 0)
  • Hapus semua pengukuran dengan kesalahan relatif tinggi ('agbd_se' / 'agbd' > 50%)
  • Hapus semua pengukuran pada lereng > 30% berdasarkan Model Ketinggian Digital (DEM) GLO-30 Copernicus

Terakhir, kita memilih semua pengukuran yang tersisa untuk jangka waktu dan wilayah yang diinginkan, lalu membuat mosaik.

var gedi = ee.ImageCollection('LARSE/GEDI/GEDI04_A_002_MONTHLY');
// Function to select the highest quality GEDI data
var qualityMask = function(image) {
  return image.updateMask(image.select('l4_quality_flag').eq(1))
      .updateMask(image.select('degrade_flag').eq(0));
};

// Function to mask unreliable GEDI measurements
// with a relative standard error > 50%
// agbd_se / agbd > 0.5
var errorMask = function(image) {
  var relative_se = image.select('agbd_se')
    .divide(image.select('agbd'));
  return image.updateMask(relative_se.lte(0.5));
};

// Function to mask GEDI measurements on slopes > 30%

var slopeMask = function(image) {
  // Use Copernicus GLO-30 DEM for calculating slope
  var glo30 = ee.ImageCollection('COPERNICUS/DEM/GLO30');

  var glo30Filtered = glo30
    .filter(ee.Filter.bounds(geometry))
    .select('DEM');

  // Extract the projection
  var demProj = glo30Filtered.first().select(0).projection();

  // The dataset consists of individual images
  // Create a mosaic and set the projection
  var elevation = glo30Filtered.mosaic().rename('dem')
    .setDefaultProjection(demProj);

  // Compute the slope
  var slope = ee.Terrain.slope(elevation);

  return image.updateMask(slope.lt(30));
};

var gediFiltered = gedi
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

var gediProjection = ee.Image(gediFiltered.first())
  .select('agbd').projection();

var gediProcessed = gediFiltered
  .map(qualityMask)
  .map(errorMask)
  .map(slopeMask);

var gediMosaic = gediProcessed.mosaic()
  .select('agbd').setDefaultProjection(gediProjection);

// Visualize the GEDI Mosaic
var gediVis = {
  min: 0,
  max: 200,
  palette: ['#edf8fb','#b2e2e2','#66c2a4','#2ca25f','#006d2c'],
  bands: ['agbd']
};

Map.addLayer(gediMosaic, gediVis, 'GEDI L4A (Filtered)', false);

Mengambil sampel ulang dan menggabungkan input

Sebelum mengambil sampel piksel untuk melatih model regresi, kita mengambil sampel ulang dan memproyeksikan ulang input ke petak piksel yang sama. Pengukuran GEDI memiliki akurasi horizontal +/- 9 m. Hal ini menimbulkan masalah saat mencocokkan nilai AGB GEDI dengan piksel Sematan Satelit. Untuk mengatasinya, kami melakukan pengambilan sampel ulang dan menggabungkan semua gambar input ke petak piksel yang lebih besar dengan nilai rata-rata dari piksel asli. Hal ini juga membantu menghilangkan derau dari data dan membantu membangun model machine learning yang lebih baik.

// Choose the grid size and projection
var gridScale = 100;
var gridProjection = ee.Projection('EPSG:3857')
  .atScale(gridScale);

// Create a stacked image with predictor and predicted variables
var stacked = embeddingsImage.addBands(gediMosaic);

//  Set the resampling mode
var stacked = stacked.resample('bilinear');

// Aggregate pixels with 'mean' statistics
var stackedResampled = stacked
  .reduceResolution({
    reducer: ee.Reducer.mean(),
    maxPixels: 1024
  })
  .reproject({
    crs: gridProjection
});

// As larger GEDI pixels contain masked original
// pixels, it has a transparency mask.
// We update the mask to remove the transparency
var stackedResampled = stackedResampled
  .updateMask(stackedResampled.mask().gt(0));

Memproyeksikan ulang dan menggabungkan piksel adalah operasi yang mahal, dan merupakan praktik yang baik untuk mengekspor gambar bertumpuk yang dihasilkan sebagai Aset dan menggunakan gambar yang telah dihitung sebelumnya dalam langkah-langkah berikutnya. Hal ini akan membantu mengatasi error waktu komputasi habis atau memori pengguna terlampaui saat bekerja dengan wilayah yang besar.

// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';
var mosaicExportImage = 'gedi_mosaic';
var mosaicExportImagePath = exportFolder + mosaicExportImage;
Export.image.toAsset({
  image: stackedResampled.clip(geometry),
  description: 'GEDI_Mosaic_Export',
  assetId: mosaicExportImagePath,
  region: geometry,
  scale: gridScale,
  maxPixels: 1e10
});

Mulai tugas ekspor dan tunggu hingga selesai. Setelah selesai, kita mengimpor Aset dan melanjutkan pembuatan model.

// Use the exported asset
var stackedResampled = ee.Image(mosaicExportImagePath);

Ekstrak fitur pelatihan

Kita telah menyiapkan data input untuk mengekstrak fitur pelatihan. Kami menggunakan band Sematan Satelit sebagai variabel dependen (prediktor) dan nilai AGBD GEDI sebagai Variabel Independen (diprediksi) dalam model regresi. Kita dapat mengekstrak nilai yang bertepatan di setiap piksel dan menyiapkan set data pelatihan. Gambar GEDI kami sebagian besar tertutup dan hanya berisi nilai pada sebagian kecil piksel. Jika kita menggunakan sample(), sebagian besar nilai yang ditampilkan akan kosong. Untuk mengatasinya, kita membuat band kelas dari mask GEDI dan menggunakan stratifiedSample() untuk memastikan kita mengambil sampel dari piksel yang tidak di-mask.

var predictors = embeddingsImage.bandNames();
var predicted = gediMosaic.bandNames().get(0);
print('predictors', predictors);
print('predicted', predicted);

var predictorImage = stackedResampled.select(predictors);
var predictedImage = stackedResampled.select([predicted]);

var classMask = predictedImage.mask().toInt().rename('class');

var numSamples = 1000;

// We set classPoints to [0, numSamples]
// This will give us 0 points for class 0 (masked areas)
// and numSample points for class 1 (non-masked areas)
var training = stackedResampled.addBands(classMask)
  .stratifiedSample({
    numPoints: numSamples,
    classBand: 'class',
    region: geometry,
    scale: gridScale,
    classValues: [0, 1],
    classPoints: [0, numSamples],
    dropNulls: true,
    tileScale: 16,
});

print('Number of Features Extracted', training.size());
print('Sample Training Feature', training.first());

Melatih model regresi

Sekarang kita siap melatih model. Banyak pengklasifikasi di Earth Engine dapat digunakan untuk tugas klasifikasi dan regresi. Karena kita ingin memprediksi nilai numerik (bukan class) – kita dapat menyetel pengklasifikasi untuk berjalan dalam mode REGRESSION dan melatih menggunakan data pelatihan. Setelah model dilatih, kita dapat membandingkan prediksi model dengan nilai input dan menghitung error kuadrat rata-rata (rmse) dan koefisien korelasi r^2 untuk memeriksa performa model.

// Use the RandomForest classifier and set the
// output mode to REGRESSION
var model = ee.Classifier.smileRandomForest(50)
  .setOutputMode('REGRESSION')
  .train({
    features: training,
    classProperty: predicted,
    inputProperties: predictors
  });

// Get model's predictions for training samples
var predicted = training.classify({
  classifier: model,
  outputName: 'agbd_predicted'
});

// Calculate RMSE
var calculateRmse = function(input) {
    var observed = ee.Array(
      input.aggregate_array('agbd'));
    var predicted = ee.Array(
      input.aggregate_array('agbd_predicted'));
    var rmse = observed.subtract(predicted).pow(2)
      .reduce('mean', [0]).sqrt().get([0]);
    return rmse;
};
var rmse = calculateRmse(predicted);
print('RMSE', rmse);

// Create a plot of observed vs. predicted values
var chart = ui.Chart.feature.byFeature({
  features: predicted.select(['agbd', 'agbd_predicted']),
  xProperty: 'agbd',
  yProperties: ['agbd_predicted'],
}).setChartType('ScatterChart')
  .setOptions({
    title: 'Aboveground Biomass Density (Mg/Ha)',
    dataOpacity: 0.8,
    hAxis: {'title': 'Observed'},
    vAxis: {'title': 'Predicted'},
    legend: {position: 'right'},
    series: {
      0: {
        visibleInLegend: false,
        color: '#525252',
        pointSize: 3,
        pointShape: 'triangle',
      },
    },
    trendlines: {
      0: {
        type: 'linear',
        color: 'black',
        lineWidth: 1,
        pointSize: 0,
        labelInLegend: 'Linear Fit',
        visibleInLegend: true,
        showR2: true
      }
    },
    chartArea: {left: 100, bottom: 100, width: '50%'},

});
print(chart);


Gambar: Nilai AGBD yang diamati vs. nilai AGBD yang diprediksi model

Membuat prediksi untuk nilai yang tidak diketahui

Setelah puas dengan modelnya, kita dapat menggunakan model terlatih untuk membuat prediksi di lokasi yang tidak diketahui dari gambar yang berisi rentang prediktor.

// We set the band name of the output image as 'agbd'
var predictedImage = stackedResampled.classify({
  classifier: model,
  outputName: 'agbd'
});

Gambar yang berisi nilai AGBD yang diprediksi di setiap piksel kini siap diekspor. Kita akan menggunakannya di bagian berikutnya untuk memvisualisasikan hasilnya.

// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';
var predictedExportImage = 'predicted_agbd';
var predictedExportImagePath = exportFolder + predictedExportImage;

Export.image.toAsset({
  image: predictedImage.clip(geometry),
  description: 'Predicted_Image_Export',
  assetId: predictedExportImagePath,
  region: geometry,
  scale: gridScale,
  maxPixels: 1e10
});

Mulai tugas ekspor dan tunggu hingga selesai. Setelah selesai, kita mengimpor Aset dan memvisualisasikan hasilnya.

var predictedImage = ee.Image(predictedExportImagePath);

// Visualize the image
var gediVis = {
  min: 0,
  max: 200,
  palette: ['#edf8fb','#b2e2e2','#66c2a4','#2ca25f','#006d2c'],
  bands: ['agbd']
};

Map.addLayer(predictedImage, gediVis, 'Predicted AGBD');


Gambar: AGBD yang diprediksi. Warna hijau yang lebih gelap menunjukkan kepadatan biomassa yang diprediksi lebih besar

Memperkirakan total biomassa

Sekarang kita memiliki nilai AGBD yang diprediksi untuk setiap piksel gambar dan nilai tersebut dapat digunakan untuk memperkirakan total stok biomassa di atas permukaan tanah (AGB) di wilayah tersebut. Namun, kita harus menghapus semua piksel yang termasuk area non-vegetasi terlebih dahulu. Kita dapat menggunakan dataset tutupan lahan ESA WorldCover dan memilih piksel yang ditumbuhi vegetasi.

// GEDI data is processed only for certain landcovers
// from Plant Functional Types (PFT) classification
// https://doi.org/10.1029/2022EA002516

// Here we use ESA WorldCover v200 product to
// select landcovers representing vegetated areas
var worldcover = ee.ImageCollection('ESA/WorldCover/v200').first();

// Aggregate pixels to the same grid as other dataset
// with 'mode' value.
// i.e. The landcover with highest occurrence within the grid
var worldcoverResampled = worldcover
  .reduceResolution({
    reducer: ee.Reducer.mode(),
    maxPixels: 1024
  })
  .reproject({
    crs: gridProjection
});

// Select grids for the following classes
// | Class Name | Value |
// | Forests    | 10    |
// | Shrubland  | 20    |
// | Grassland  | 30    |
// | Cropland   | 40    |
// | Mangroves  | 95    |
var landCoverMask = worldcoverResampled.eq(10)
    .or(worldcoverResampled.eq(20))
    .or(worldcoverResampled.eq(30))
    .or(worldcoverResampled.eq(40))
    .or(worldcoverResampled.eq(95));

var predictedImageMasked = predictedImage
  .updateMask(landCoverMask);
Map.addLayer(predictedImageMasked, gediVis, 'Predicted AGBD (Masked)');

Satuan nilai AGBD GEDI adalah megagram per hektar (Mg/ha). Untuk mendapatkan total AGB, kita mengalikan setiap piksel dengan luasnya dalam hektar dan menjumlahkan nilainya.

var pixelAreaHa = ee.Image.pixelArea().divide(10000);
var predictedAgb = predictedImageMasked.multiply(pixelAreaHa);

var stats = predictedAgb.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: geometry,
  scale: gridScale,
  maxPixels: 1e10,
  tileScale: 16
});

// Result is a dictionary with key for each band
var totalAgb = stats.getNumber('agbd');

print('Total AGB (Mg)', totalAgb);

Coba skrip lengkap untuk tutorial ini di Editor Kode Earth Engine.