Regresión con el conjunto de datos de incorporación de satélite

Editar en GitHub

Informa un problema

Historial de la página

Autor(es): spatialthoughts

Los instructivos que aporta la comunidad de desarrolladores de Earth Engine no forman parte de la documentación oficial del producto de Earth Engine.

Este instructivo forma parte de una serie de instructivos sobre el conjunto de datos de Satellite Embedding. Consulta también Introducción, Clasificación no supervisada, Clasificación supervisada y Búsqueda de similitud.

Los campos de incorporación se pueden usar como predictores o entradas de atributos para la regresión de la misma manera en que se usan para la clasificación.

En este instructivo, aprenderemos a usar las capas de campos de incorporación de 64 dimensiones como entradas para un análisis de regresión múltiple que predice la biomasa sobre el suelo (AGB).

La misión Global Ecosystem Dynamics Investigation (GEDI) de la NASA recopila mediciones de LIDAR a lo largo de transectos terrestres con una resolución espacial de 30 m en intervalos de 60 m. Usaremos el conjunto de datos GEDI L4A Raster Aboveground Biomass Density, que contiene estimaciones puntuales de la densidad de biomasa sobre el suelo (AGBD) que se utilizarán como variable predictiva en el modelo de regresión.

Selecciona una región

Comencemos por definir una región de interés. En este instructivo, elegiremos una región en los Ghats occidentales de la India y definiremos un polígono como la variable de geometría. Como alternativa, puedes usar las herramientas de dibujo en el editor de código para dibujar un polígono alrededor de la región de interés que se guardará como la variable de geometría en las importaciones. También usamos el mapa base de satélite, que facilita la ubicación de las áreas con vegetación.

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: Selección del área de interés para la predicción de biomasa sobre el suelo

Selecciona un período

Elige un año para el que queremos ejecutar la regresión. Recuerda que las incorporaciones de satélite se agregan en intervalos anuales, por lo que definimos el período para todo el año.

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

Prepara el conjunto de datos de Satellite Embedding

Las imágenes de Satellite Embedding de 64 bandas se usarán como predictor para la regresión. Cargamos el conjunto de datos de Satellite Embedding y filtramos las imágenes del año y la región elegidos.

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

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

Las imágenes de Satellite Embedding se dividen en cuadrículas y se publican en la proyección de las zonas UTM de la tarjeta. Como resultado, obtenemos varias tarjetas de Satellite Embedding que cubren la región de interés. Para obtener una sola imagen, debemos crear un mosaico con ellas. En Earth Engine, a un mosaico de imágenes de entrada se le asigna la proyección predeterminada, que es WGS84 con una escala de 1 grado. Como agregaremos y reproyectaremos este mosaico más adelante en el instructivo, es útil conservar la proyección original. Podemos extraer la información de proyección de una de las tarjetas y establecerla en el mosaico con la función 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 el mosaico L4A de GEDI

Dado que las estimaciones de biomasa de GEDI se usarán para entrenar nuestro modelo de regresión, es fundamental filtrar los datos de GEDI no válidos o poco confiables antes de usarlos. Aplicamos varias máscaras para quitar las mediciones potencialmente erróneas.

• Quita todas las mediciones que no cumplan con el requisito de calidad (l4_quality_flag = 0 y degrade_flag > 0).
• Quita todas las mediciones con un error relativo alto ('agbd_se' / 'agbd' > 50%).
• Quita todas las mediciones en pendientes superiores al 30% según el Modelo digital de elevación (MDE) GLO-30 de Copernicus.

Por último, seleccionamos todas las mediciones restantes para el período y la región de interés, y creamos 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);

Cómo volver a muestrear y agregar entradas

Antes de tomar muestras de píxeles para entrenar un modelo de regresión, volvemos a muestrear y proyectar las entradas en la misma cuadrícula de píxeles. Las mediciones del GEDI tienen una precisión horizontal de +/- 9 m. Esto es problemático cuando se comparan los valores de AGB de GEDI con los píxeles de Satellite Embedding. Para superar este problema, volvemos a muestrear y agregamos todas las imágenes de entrada en una cuadrícula de píxeles más grande con valores medios de los píxeles originales. Esto también ayuda a quitar el ruido de los datos y a crear un mejor modelo de aprendizaje automático.

// 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 reproyección y la agregación de píxeles son operaciones costosas, por lo que es una práctica recomendada exportar la imagen apilada resultante como un recurso y usar la imagen precalculada en los pasos posteriores. Esto ayudará a superar los errores de tiempo de espera agotado para el cálculo o se excedió la memoria del usuario cuando se trabaja con regiones grandes.

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

Inicia la tarea de exportación y espera a que finalice. Una vez que terminemos, importaremos el activo y continuaremos creando el modelo.

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

Extrae atributos de entrenamiento

Tenemos nuestros datos de entrada listos para extraer atributos de entrenamiento. Utilizamos las bandas de Satellite Embedding como variables dependientes (predictoras) y los valores de AGBD de GEDI como variable independiente (predicha) en el modelo de regresión. Podemos extraer los valores coincidentes en cada píxel y preparar nuestro conjunto de datos de entrenamiento. Nuestra imagen de GEDI está enmascarada en su mayor parte y contiene valores solo en un pequeño subconjunto de píxeles. Si usamos sample(), se devolverán valores casi vacíos. Para superar este problema, creamos una banda de clases a partir de la máscara de GEDI y usamos stratifiedSample() para asegurarnos de que tomamos muestras de los píxeles no enmascarados.

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

Entrena un modelo de regresión

Ahora sí podemos entrenar el modelo. Muchos clasificadores en Earth Engine se pueden usar para tareas de clasificación y regresión. Dado que queremos predecir un valor numérico (en lugar de una clase), podemos configurar el clasificador para que se ejecute en el modo REGRESSION y entrenarlo con los datos de entrenamiento. Una vez que se entrena el modelo, podemos comparar su predicción con los valores de entrada y calcular el error cuadrático medio (rmse) y el coeficiente de correlación r^2 para verificar el rendimiento del modelo.

// 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: Comparación entre los valores de AGBD observados y los predichos por el modelo

Genera predicciones para valores desconocidos

Una vez que estemos satisfechos con el modelo, podemos usar el modelo entrenado para generar predicciones en ubicaciones desconocidas a partir de la imagen que contiene bandas de predictores.

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

La imagen que contiene los valores de AGBD predichos en cada píxel ya está lista para exportarse. Usaremos esto en la próxima parte para visualizar los resultados.

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

Inicia la tarea de exportación y espera a que finalice. Una vez que terminemos, importaremos el activo y visualizaremos los resultados.

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 prevista. Los tonos de verde más oscuros indican una mayor densidad de biomasa prevista

Estima la biomasa total

Ahora tenemos valores de AGBD predichos para cada píxel de la imagen, que se pueden usar para estimar el stock total de biomasa sobre el suelo (AGB) en la región. Sin embargo, primero debemos quitar todos los píxeles que pertenecen a áreas no vegetadas. Podemos usar el conjunto de datos de cobertura terrestre ESA WorldCover y seleccionar los píxeles con vegetación.

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

Las unidades de los valores de AGBD de GEDI son megagramos por hectárea (Mg/ha). Para obtener el AGB total, multiplicamos cada píxel por su área en hectáreas y sumamos sus valores.

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

Prueba la secuencia de comandos completa de este instructivo en el editor de código de Earth Engine.