Regressão com conjunto de dados de embedding de satélite

Editar no GitHub

Informar problema

Histórico da página

Autor(es): spatialthoughts

Os tutoriais contribuídos pela comunidade de desenvolvedores do Earth Engine não fazem parte da documentação oficial do produto.

Este tutorial faz parte de uma série sobre o conjunto de dados de incorporação de satélite. Consulte também Introdução, Classificação não supervisionada, Classificação supervisionada e Pesquisa de similaridade.

Os campos de embedding podem ser usados como entradas/preditores de recursos para regressão da mesma forma que são usados para classificação.

Neste tutorial, vamos aprender a usar as camadas de campo de incorporação 64D como entradas para uma análise de regressão múltipla que prevê a biomassa acima do solo (AGB, na sigla em inglês).

A missão Global Ecosystem Dynamics Investigation (GEDI) da NASA coleta medições de LIDAR ao longo de transectos terrestres com resolução espacial de 30 m em intervalos de 60 m. Vamos usar o conjunto de dados GEDI L4A Raster Aboveground Biomass Density, que contém estimativas pontuais de densidade de biomassa acima do solo (AGBD, na sigla em inglês) e será usado como a variável prevista no modelo de regressão.

Selecione uma região

Vamos começar definindo uma região de interesse. Para este tutorial, vamos escolher uma região nos Gates Ocidentais da Índia e definir um polígono como a variável de geometria. Como alternativa, use as ferramentas de desenho no editor de código para desenhar um polígono ao redor da região de interesse, que será salva como a variável de geometria nas importações. Também usamos o mapa de base de satélite, que facilita a localização de áreas com vegetação.

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: seleção da área de interesse para previsão de biomassa acima do solo

Selecione um período

Escolha um ano para executar a regressão. Lembre-se de que as incorporações de satélite são agregadas em intervalos anuais. Portanto, definimos o período para o ano inteiro.

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

Preparar o conjunto de dados de embedding de satélite

As imagens de incorporação de satélite de 64 bandas serão usadas como preditor para a regressão. Carregamos o conjunto de dados de incorporação de satélite e filtramos as imagens do ano e da região escolhidos.

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

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

As imagens de incorporação de satélite são divididas em blocos e veiculadas na projeção das zonas UTM do bloco. Como resultado, recebemos vários blocos de incorporação de satélite que cobrem a região de interesse. Para ter uma única imagem, precisamos criar um mosaico. No Earth Engine, um mosaico de imagens de entrada recebe a projeção padrão, que é WGS84 com uma escala de 1 grau. Como vamos agregar e reprojetar esse mosaico mais tarde no tutorial, é útil manter a projeção original. Podemos extrair as informações de projeção de um dos blocos e defini-las no mosaico usando a função 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);

Preparar o mosaico L4A do GEDI

Como as estimativas de biomassa do GEDI serão usadas para treinar nosso modelo de regressão, é fundamental filtrar os dados inválidos ou não confiáveis antes de usá-los. Aplicamos várias máscaras para remover medições potencialmente errôneas.

• Remova todas as medições que não atendem ao requisito de qualidade (l4_quality_flag = 0 e degrade_flag > 0)
• Remova todas as medições com erro relativo alto ("agbd_se" / "agbd" > 50%)
• Remover todas as medições em inclinações > 30% com base no Modelo Digital de Elevação (MDE) GLO-30 do Copernicus

Por fim, selecionamos todas as medições restantes para o período e a região de interesse e criamos um 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);

Reamostrar e agregar entradas

Antes de fazer a amostragem de pixels para treinar um modelo de regressão, fazemos uma nova amostragem e reprojetamos as entradas para a mesma grade de pixels. As medições do GEDI têm uma precisão horizontal de +/- 9 m. Isso é problemático ao corresponder os valores da AGB do GEDI aos pixels de incorporação de satélite. Para resolver isso, fazemos uma nova amostragem e agregamos todas as imagens de entrada em uma grade de pixels maior com valores médios dos pixels originais. Isso também ajuda a remover ruídos dos dados e a criar um modelo de aprendizado de máquina melhor.

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

A reprojeção e a agregação de pixels são operações caras. Por isso, é recomendado exportar a imagem combinada resultante como um recurso e usar a imagem pré-calculada nas etapas subsequentes. Isso ajuda a superar erros de tempo limite de computação excedido ou memória do usuário excedida ao trabalhar com regiões 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
});

Inicie a tarefa de exportação e aguarde a conclusão. Depois disso, importamos o recurso e continuamos criando o modelo.

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

Extrair atributos de treinamento

Os dados de entrada estão prontos para extrair recursos de treinamento. Usamos as bandas de embedding do satélite como variáveis dependentes (preditoras) e os valores de AGBD do GEDI como variáveis independentes (previstas) no modelo de regressão. Podemos extrair os valores coincidentes em cada pixel e preparar nosso conjunto de dados de treinamento. Nossa imagem GEDI é principalmente mascarada e contém valores apenas em um pequeno subconjunto de pixels. Se usarmos sample(), ele vai retornar principalmente valores vazios. Para resolver isso, criamos uma banda de classe com a máscara do GEDI e usamos stratifiedSample() para garantir que façamos a amostragem dos pixels não mascarados.

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

Treinar um modelo de regressão

Agora estamos prontos para treinar o modelo. Muitos classificadores no Earth Engine podem ser usados para tarefas de classificação e regressão. Como queremos prever um valor numérico (em vez de uma classe), podemos definir o classificador para ser executado no modo REGRESSION e treinar usando os dados de treinamento. Depois que o modelo é treinado, podemos comparar a previsão dele com os valores de entrada e calcular o erro quadrático médio (rmse) e o coeficiente de correlação r^2 para verificar o desempenho do 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: valores de AGBD observados x previstos pelo modelo

Gerar previsões para valores desconhecidos

Quando estivermos satisfeitos com o modelo, poderemos usá-lo para gerar previsões em locais desconhecidos na imagem com intervalos de previsão.

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

A imagem com os valores previstos de AGBD em cada pixel está pronta para exportação. Vamos usar isso na próxima parte para visualizar os 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
});

Inicie a tarefa de exportação e aguarde a conclusão. Depois disso, importamos o recurso e visualizamos os 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 previsto. Os tons mais escuros indicam uma densidade de biomassa prevista maior.

Estimar a biomassa total

Agora temos valores previstos de AGBD para cada pixel da imagem, que podem ser usados para estimar o estoque total de biomassa acima do solo (AGB, na sigla em inglês) na região. Mas primeiro precisamos remover todos os pixels pertencentes a áreas não vegetadas. Podemos usar o conjunto de dados de cobertura da terra ESA WorldCover e selecionar pixels vegetados.

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

As unidades dos valores de AGBD do GEDI são megagramas por hectare (Mg/ha). Para saber o AGB total, multiplicamos cada pixel pela área em hectares e somamos os 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);

Teste o script completo deste tutorial no editor de código do Earth Engine.