Regressione con il set di dati di incorporamento satellitare

Modifica su GitHub
Segnala problema
Cronologia della pagina
Questo tutorial fa parte di una serie di tutorial sul set di dati Satellite Embedding. Vedi anche Introduzione, Classificazione non supervisionata, Classificazione supervisionata e Ricerca per somiglianza.

I campi di incorporamento possono essere utilizzati come input/predittori delle funzionalità per la regressione nello stesso modo in cui vengono utilizzati per la classificazione.

In questo tutorial impareremo a utilizzare i livelli del campo di incorporamento 64D come input per un'analisi di regressione multipla che prevede la biomassa fuori terra (AGB).

La missione Global Ecosystem Dynamics Investigation (GEDI) della NASA raccoglie misurazioni LIDAR lungo i transetti a terra con una risoluzione spaziale di 30 m a intervalli di 60 m. Utilizzeremo il set di dati GEDI L4A Raster Aboveground Biomass Density contenente stime puntuali della densità di biomassa epigea (AGBD), che verrà utilizzato come variabile prevista nel modello di regressione.

Seleziona una regione

Iniziamo definendo una regione di interesse. Per questo tutorial, sceglieremo una regione nei Ghati occidentali dell'India e definiremo un poligono come variabile di geometria. In alternativa, puoi utilizzare gli strumenti di disegno nell'editor di codice per disegnare un poligono intorno alla regione di interesse che verrà salvato come variabile di geometria nelle importazioni. Utilizziamo anche la mappa di base satellitare, che consente di individuare facilmente le aree con vegetazione.

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


Figura: selezione dell'area di interesse per la previsione della biomassa fuori terra

Seleziona un periodo di tempo

Scegli un anno per cui eseguire la regressione. Tieni presente che gli incorporamenti satellitari vengono aggregati a intervalli annuali, quindi definiamo il periodo per l'intero anno.

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

Prepara il set di dati di incorporamento satellitare

Le immagini di incorporamento satellitare a 64 bande verranno utilizzate come predittore per la regressione. Carichiamo il set di dati Satellite Embedding, filtriamo le immagini per l'anno e la regione scelti.

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

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

Le immagini incorporate del satellite vengono suddivise in riquadri e visualizzate nella proiezione per le zone UTM del riquadro. Di conseguenza, otteniamo più riquadri di incorporamento satellitare che coprono la regione di interesse. Per ottenere un'unica immagine, dobbiamo creare un mosaico. In Earth Engine, a un mosaico di immagini di input viene assegnata la proiezione predefinita, ovvero WGS84 con una scala di 1 grado. Poiché questo mosaico verrà aggregato e riproiettato più avanti nel tutorial, è utile conservare la proiezione originale. Possiamo estrarre le informazioni sulla proiezione da uno dei riquadri e impostarle sul mosaico utilizzando la funzione 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);

Prepara il mosaico GEDI L4A

Poiché le stime della biomassa GEDI verranno utilizzate per addestrare il nostro modello di regressione, è fondamentale filtrare i dati GEDI non validi o inaffidabili prima di utilizzarli. Applichiamo diverse maschere per rimuovere le misurazioni potenzialmente errate.

  • Rimuovi tutte le misurazioni che non soddisfano il requisito di qualità (l4_quality_flag = 0 e degrade_flag > 0)
  • Rimuovi tutte le misurazioni con un errore relativo elevato ("agbd_se" / "agbd" > 50%)
  • Rimuovi tutte le misurazioni su pendenze > 30% in base al modello digitale di elevazione (DEM) Copernicus GLO-30

Infine, selezioniamo tutte le misurazioni rimanenti per il periodo di tempo e la regione di interesse e creiamo un mosaico.

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

Esegui il ricampionamento e l'aggregazione degli input

Prima di campionare i pixel per addestrare un modello di regressione, ricampioniamo e riproiettiamo gli input nella stessa griglia di pixel. Le misurazioni GEDI hanno una precisione orizzontale di +/- 9 m. Questo è problematico quando si abbinano i valori GEDI AGB ai pixel di incorporamento satellitare. Per superare questo problema, eseguiamo il ricampionamento e l'aggregazione di tutte le immagini di input in una griglia di pixel più grande con i valori medi dei pixel originali. Ciò contribuisce anche a eliminare il rumore dai dati e a creare un modello di machine learning migliore.

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

La riproiezione e l'aggregazione dei pixel sono operazioni costose ed è una buona pratica esportare l'immagine in pila risultante come asset e utilizzare l'immagine pre-elaborata nei passaggi successivi. In questo modo, potrai superare gli errori Calcolo scaduto o Memoria utente superata quando lavori con regioni di grandi dimensioni.

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

Avvia l'attività di esportazione e attendi il completamento. Al termine, importiamo l'asset e continuiamo a creare il modello.

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

Estrai le caratteristiche di addestramento

Abbiamo i dati di input pronti per l'estrazione delle caratteristiche di addestramento. Utilizziamo le bande di incorporamento del satellite come variabili dipendenti (predittori) e i valori GEDI AGBD come variabile indipendente (prevista) nel modello di regressione. Possiamo estrarre i valori coincidenti in ogni pixel e preparare il nostro set di dati di addestramento. La nostra immagine GEDI è per lo più mascherata e contiene valori solo in un piccolo sottoinsieme di pixel. Se utilizziamo sample(), verranno restituiti per lo più valori vuoti. Per superare questo problema, creiamo una banda di classe dalla maschera GEDI e utilizziamo stratifiedSample() per assicurarci di campionare i pixel non mascherati.

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

Addestra un modello di regressione

Ora possiamo addestrare il modello. Molti classificatori in Earth Engine possono essere utilizzati per attività di classificazione e regressione. Poiché vogliamo prevedere un valore numerico (anziché una classe), possiamo impostare il classificatore per l'esecuzione in modalità REGRESSION e addestrarlo utilizzando i dati di addestramento. Una volta addestrato il modello, possiamo confrontare la previsione del modello con i valori di input e calcolare l'errore quadratico medio (rmse) e il coefficiente di correlazione r^2 per verificare le prestazioni del modello.

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


Figura: valori AGBD osservati e previsti dal modello

Generare previsioni per valori sconosciuti

Una volta che il modello ci soddisfa, possiamo utilizzarlo per generare previsioni in località sconosciute dall'immagine contenente le bande del predittore.

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

L'immagine contenente i valori AGBD previsti per ogni pixel è ora pronta per l'esportazione. Lo utilizzeremo nella parte successiva per visualizzare i risultati.

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

Avvia l'attività di esportazione e attendi il completamento. Al termine, importiamo l'asset e visualizziamo i risultati.

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


Figura: AGBD previsto. I verdi più scuri indicano una maggiore densità di biomassa prevista

Stima della biomassa totale

Ora abbiamo valori AGBD previsti per ogni pixel dell'immagine, che possono essere utilizzati per stimare lo stock totale di biomassa epigea (AGB) nella regione. Ma prima dobbiamo rimuovere tutti i pixel appartenenti ad aree non vegetate. Possiamo utilizzare il set di dati di copertura del suolo ESA WorldCover e selezionare i pixel vegetati.

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

Le unità dei valori GEDI AGBD sono megagrammi per ettaro (Mg/ha). Per ottenere l'AGB totale, moltiplichiamo ogni pixel per la sua area in ettari e sommiamo i valori.

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

Prova lo script completo per questo tutorial nell'editor di codice di Earth Engine.