Uydu Yerleştirme Veri Kümesi ile Regresyon

GitHub'da düzenle
Sorun bildir
Sayfa geçmişi
Bu eğitim, Satellite Embedding veri kümesiyle ilgili bir dizi eğitimin parçasıdır. Ayrıca Giriş, Gözetimsiz Sınıflandırma, Gözetimli Sınıflandırma ve Benzerlik Arama eğitimlerine de göz atın.

Yerleştirme alanları, sınıflandırmada kullanıldığı şekilde regresyon için özellik girişi/tahminci olarak kullanılabilir.

Bu eğitimde, yer üstü biyokütlesini (AGB) tahmin eden çoklu regresyon analizine giriş olarak 64D gömme alan katmanlarının nasıl kullanılacağını öğreneceğiz.

NASA'nın Global Ecosystem Dynamics Investigation (GEDI) misyonu, 60 m aralıklarla 30 m mekansal çözünürlükte yer transekleri boyunca LIDAR ölçümleri toplar. Regresyon modelinde tahmin edilen değişken olarak kullanılacak, yer üstü biyokütle yoğunluğunun (AGBD) nokta tahminlerini içeren GEDI L4A Raster Aboveground Biomass Density veri kümesini kullanacağız.

Bölge seçin

İşe ilgi alanını tanımlayarak başlayalım. Bu eğitimde, Hindistan'ın Batı Ghats bölgesinde bir bölge seçecek ve geometri değişkeni olarak bir poligon tanımlayacağız. Alternatif olarak, ilgi alanının etrafına bir poligon çizmek için kod düzenleyicideki Çizim Araçları'nı kullanabilirsiniz. Bu poligon, içe aktarma işlemlerinde geometri değişkeni olarak kaydedilir. Ayrıca, bitki örtüsü olan alanları bulmayı kolaylaştıran uydu temel haritasını da kullanırız.

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


Şekil: Yer üstü biyokütle tahmini için ilgi alanını seçme

Bir zaman aralığı seçin

Regresyonu çalıştırmak istediğimiz yılı seçin. Uydu yerleştirmelerinin yıllık aralıklarla toplandığını unutmayın. Bu nedenle, dönemi tüm yıl olarak tanımlarız.

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

Uydu yerleştirme veri kümesini hazırlama

64 bantlı uydu yerleştirme görüntüleri, regresyon için tahmin edici olarak kullanılır. Uydu Yerleştirme veri kümesini yüklüyoruz, seçilen yıl ve bölgeye ait görüntüleri filtreliyoruz.

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

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

Uydu yerleştirme resimleri, karolar halinde ızgaralara bölünür ve karo için UTM bölgelerinin projeksiyonunda sunulur. Bu nedenle, ilgilenilen bölgeyi kapsayan birden fazla uydu yerleştirme kutucuğu elde ederiz. Tek bir görüntü elde etmek için bunları mozaik haline getirmemiz gerekir. Earth Engine'de, giriş resimlerinin mozaik görüntüsüne varsayılan projeksiyon atanır. Bu projeksiyon, 1 derecelik ölçeğe sahip WGS84'tür. Bu mozaik daha sonra eğitimde toplanıp yeniden yansıtılacağından orijinal yansıtmayı korumak faydalıdır. Yansıtma bilgilerini kutulardan birinden çıkarabilir ve setDefaultProjection() işlevini kullanarak mozaik üzerinde ayarlayabiliriz.

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

GEDI L4A mozaikini hazırlama

GEDI biyokütle tahminleri regresyon modelimizi eğitmek için kullanılacağından, kullanmadan önce geçersiz veya güvenilir olmayan GEDI verilerini filtrelemek çok önemlidir. Hatalı olabilecek ölçümleri kaldırmak için çeşitli maskeler uygularız.

  • Kalite şartını karşılamayan tüm ölçümleri kaldırın (l4_quality_flag = 0 ve degrade_flag > 0)
  • Yüksek bağıl hataya sahip tüm ölçümleri kaldırın ("agbd_se" / "agbd" > %50).
  • Copernicus GLO-30 Dijital Yükseklik Modu'na (DEM) göre% 30'dan yüksek eğimlerdeki tüm ölçümleri kaldırın.

Son olarak, ilgilenilen dönem ve bölge için kalan tüm ölçümleri seçip mozaik oluştururuz.

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

Girişleri yeniden örnekleme ve toplama

Bir regresyon modelini eğitmek için pikselleri örneklemeden önce girişleri aynı piksel ızgarasına yeniden örnekleyip yeniden yansıtıyoruz. GEDI ölçümleri +/- 9 m yatay doğruluğa sahiptir. Bu durum, GEDI AGB değerleri ile Satellite Embedding pikselleri eşleştirilirken sorunlara yol açar. Bunun üstesinden gelmek için tüm giriş resimlerini yeniden örnekleyip orijinal piksellerdeki ortalama değerlerle daha büyük bir piksel ızgarasında topluyoruz. Bu, verilerdeki gürültüyü kaldırmaya ve daha iyi bir makine öğrenimi modeli oluşturmaya da yardımcı olur.

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

Piksellerin yeniden yansıtılması ve toplanması maliyetli bir işlemdir. Bu nedenle, elde edilen yığılmış görüntüyü öğe olarak dışa aktarmak ve sonraki adımlarda önceden hesaplanmış görüntüyü kullanmak iyi bir uygulamadır. Bu, büyük bölgelerle çalışırken işlem zaman aşımına uğradı veya kullanıcı belleği aşıldı hatalarının giderilmesine yardımcı olur.

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

Dışa aktarma görevini başlatın ve tamamlanmasını bekleyin. İşlem tamamlandıktan sonra öğeyi içe aktarıp modeli oluşturmaya devam ederiz.

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

Eğitim özelliklerini çıkarma

Eğitim özelliklerini çıkarmak için giriş verilerimiz hazır. Regresyon modelinde bağımlı değişkenler (tahmin ediciler) olarak uydu yerleştirme bantlarını, bağımsız değişken (tahmin edilen) olarak ise GEDI AGBD değerlerini kullanırız. Her pikseldeki çakışan değerleri çıkarıp eğitim veri kümemizi hazırlayabiliriz. GEDI görüntümüzün çoğu maskelenmiş durumda ve yalnızca küçük bir piksel alt kümesinde değerler içeriyor. sample() kullanırsak çoğunlukla boş değerler döndürülür. Bu sorunu aşmak için GEDI maskesinden bir sınıf bandı oluştururuz ve maskelenmemiş piksellerden örnekleme yaptığımızdan emin olmak için stratifiedSample() kullanırız.

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

Regresyon modeli eğitme

Artık modeli eğitmeye hazırız. Earth Engine'deki birçok sınıflandırıcı, hem sınıflandırma hem de regresyon görevleri için kullanılabilir. Bir sınıf yerine sayısal bir değer tahmin etmek istediğimiz için sınıflandırıcıyı REGRESSION modunda çalışacak şekilde ayarlayabilir ve eğitim verilerini kullanarak eğitebiliriz. Model eğitildikten sonra, modelin performansını kontrol etmek için modelin tahmini ile giriş değerlerini karşılaştırabilir ve ortalama kare hatayı (rmse) ve korelasyon katsayısını (r^2) hesaplayabiliriz.

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


Şekil: Gözlemlenen ve modelin tahmin ettiği AGBD değerleri

Bilinmeyen değerler için tahmin oluşturma

Modelden memnun kaldığımızda, eğitilmiş modeli kullanarak tahmin bantlarını içeren resimdeki bilinmeyen konumlarda tahminler oluşturabiliriz.

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

Her pikselde tahmin edilen AGBD değerlerini içeren resim artık dışa aktarılmaya hazır. Sonuçları görselleştirmek için bunu bir sonraki bölümde kullanacağız.

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

Dışa aktarma görevini başlatın ve tamamlanmasını bekleyin. İşlem tamamlandıktan sonra öğeyi içe aktarıp sonuçları görselleştiririz.

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


Şekil: Tahmini AGBD. Daha koyu yeşiller, daha yüksek tahmini biyokütle yoğunluğunu gösterir.

Toplam biyokütleyi tahmin etme

Artık resmin her pikseli için tahmin edilen AGBD değerlerimiz var. Bu değerler, bölgedeki toplam yer üstü biyokütle (AGB) stokunu tahmin etmek için kullanılabilir. Ancak önce bitki örtüsü olmayan alanlara ait tüm pikselleri kaldırmamız gerekir. ESA WorldCover arazi örtüsü veri kümesini kullanabilir ve bitki örtüsü olan pikselleri seçebiliriz.

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

GEDI AGBD değerlerinin birimleri hektar başına megagramdır (Mg/ha). Toplam AGB'yi elde etmek için her pikseli hektar cinsinden alanıyla çarpar ve değerlerini toplarız.

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

Bu eğitimin tam komut dosyasını Earth Engine Kod Düzenleyici'de deneyin.