Clasificación no supervisada con el conjunto de datos de Satellite Embedding

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 supervisada, Regresión y Búsqueda de similitud.

En el instructivo anterior (Introducción), vimos cómo los embeddings de satélite capturan las trayectorias anuales de las observaciones satelitales y las variables climáticas. Esto hace que el conjunto de datos sea fácilmente utilizable para mapear cultivos sin necesidad de modelar la fenología de los cultivos. La asignación de tipos de cultivos es una tarea desafiante que, por lo general, requiere modelar la fenología de los cultivos y recopilar muestras de campo de todos los cultivos que se producen en la región.

En este instructivo, adoptaremos un enfoque de clasificación no supervisada para la asignación de cultivos que nos permitirá realizar esta tarea compleja sin depender de etiquetas de campo. Este método aprovecha el conocimiento local de la región junto con las estadísticas agregadas de los cultivos, que están disponibles para muchas partes del mundo.

Selecciona una región

En este instructivo, crearemos un mapa de tipos de cultivos para el condado de Cerro Gordo, Iowa. Este condado se encuentra en el cinturón de maíz de Estados Unidos, que tiene dos cultivos principales: maíz y soja. Este conocimiento local es importante y nos ayudará a decidir los parámetros clave de nuestro modelo.

Comencemos por obtener el límite del condado elegido.

// Select the region
// Cerro Gordo County, Iowa
var counties = ee.FeatureCollection('TIGER/2018/Counties');

var selected = counties
  .filter(ee.Filter.eq('GEOID', '19033'));
var geometry = selected.geometry();
Map.centerObject(geometry, 12);
Map.addLayer(geometry, {color: 'red'}, 'Selected Region', false);

Figura: Región seleccionada

Prepara el conjunto de datos de Satellite Embedding

A continuación, cargamos el conjunto de datos de Satellite Embedding, filtramos las imágenes del año elegido y creamos un mosaico.

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

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

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

var embeddingsImage = filteredembeddings.mosaic();

Crea una máscara de tierras de cultivo

Para nuestro modelado, debemos excluir las áreas que no son de cultivo. Existen muchos conjuntos de datos globales y regionales que se pueden usar para crear una máscara de recorte. ESA WorldCover o GFSAD Global Cropland Extent Product son buenas opciones para los conjuntos de datos de tierras de cultivo globales. Una incorporación más reciente es el producto ESA WorldCereal Active Cropland, que incluye el mapeo estacional de las tierras de cultivo activas. Dado que nuestra región se encuentra en EE.UU., podemos usar un conjunto de datos regional más preciso, USDA NASS Cropland Data Layers (CDL), para obtener una máscara de cultivo.

// Use Cropland Data Layers (CDL) to obtain cultivated cropland
var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var croplandMask = cdl.select('cultivated').eq(2).rename('cropmask');

// Visualize the crop mask
var croplandMaskVis = {min: 0, max: 1, palette: ['white', 'green']};
Map.addLayer(croplandMask.clip(geometry), croplandMaskVis, 'Crop Mask');

Figura: Región seleccionada con máscara de tierra de cultivo

Extrae muestras de entrenamiento

Aplicamos la máscara de recorte a la composición de mosaicos de la incorporación. Ahora, nos quedan todos los píxeles que representan las tierras de cultivo del condado.

// Mask all non-cropland pixels
var clusterImage = embeddingsImage.updateMask(croplandMask);

Debemos tomar la imagen de Satellite Embedding y obtener muestras aleatorias para entrenar un modelo de clustering. Dado que nuestra región de interés contiene muchos píxeles enmascarados, un muestreo aleatorio simple puede generar muestras con valores nulos. Para asegurarnos de que podemos extraer la cantidad deseada de muestras no nulas, usamos el muestreo estratificado para obtener la cantidad deseada de muestras en las áreas sin máscara.

// Stratified random sampling
var training = clusterImage.addBands(croplandMask).stratifiedSample({
  numPoints: 1000,
  classBand: 'cropmask',
  region: geometry,
  scale: 10,
  tileScale: 16,
  seed: 100,
  dropNulls: true,
  geometries: true
});

Exporta la muestra a un activo (opcional)

Nota: Este ejemplo está diseñado para funcionar de forma interactiva, pero es posible que no se escale a regiones más grandes o a una mayor cantidad de muestras. En este caso, puedes usar una exportación para solucionar los problemas de escalamiento.

La extracción de muestras es una operación costosa en términos de procesamiento, por lo que es una práctica recomendada exportar las muestras de entrenamiento extraídas como recursos y usar los recursos exportados en los pasos posteriores. Esto ayudará a superar los errores de tiempo de espera agotado para el cálculo o memoria del usuario excedida cuando se trabaja con regiones grandes.

Inicia la tarea de exportación y espera a que finalice antes de continuar.

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

var samplesExportFc = 'cluster_training_samples';
var samplesExportFcPath = exportFolder + samplesExportFc;

Export.table.toAsset({
  collection: training,
  description: 'Cluster_Training_Samples',
  assetId: samplesExportFcPath
});

Una vez que se complete la tarea de exportación, podremos leer las muestras extraídas en nuestro código como una colección de atributos.

// Use the exported asset
var training = ee.FeatureCollection(samplesExportFcPath);

Visualiza las muestras

Ya sea que hayas ejecutado tu muestra de forma interactiva o la hayas exportado a una colección de entidades, ahora tendrás una variable de entrenamiento con tus puntos de muestra. Imprimamos la primera muestra para inspeccionar y agregar nuestros puntos de entrenamiento a Map.

print('Extracted sample', training.first());
Map.addLayer(training, {color: 'blue'}, 'Extracted Samples', false);

Figura: Muestras aleatorias extraídas para la agrupación en clústeres

Realiza el agrupamiento en clústeres no supervisado

Ahora podemos entrenar un agrupador en clústeres y agrupar los vectores de embedding de 64 dimensiones en una cantidad elegida de clústeres distintos. Según nuestro conocimiento local, existen dos tipos principales de cultivos que representan la mayoría del área, y varios otros cultivos cubren la fracción restante. Podemos realizar una agrupación no supervisada en la incorporación de satélite para obtener clústeres de píxeles que tengan patrones y trayectorias temporales similares. Los píxeles con características espectrales y espaciales similares, junto con una fenología similar, se agruparán en el mismo clúster.

El ee.Clusterer.wekaCascadeKMeans() nos permite especificar una cantidad mínima y máxima de clústeres, y encontrar la cantidad óptima de clústeres según los datos de entrenamiento. Aquí es donde nuestro conocimiento local resulta útil para decidir la cantidad mínima y máxima de clústeres. Dado que esperamos algunos tipos distintos de clústeres (maíz, soja y varios otros cultivos), podemos usar 4 como la cantidad mínima de clústeres y 5 como la cantidad máxima de clústeres. Es posible que debas experimentar con estas cifras para descubrir qué funciona mejor en tu región.

// Cluster the Satellite Embedding Image
var minClusters = 4;
var maxClusters = 5;

var clusterer = ee.Clusterer.wekaCascadeKMeans({
  minClusters: minClusters, maxClusters: maxClusters}).train({
  features: training,
  inputProperties: clusterImage.bandNames()
});

var clustered = clusterImage.cluster(clusterer);
Map.addLayer(clustered.randomVisualizer().clip(geometry), {}, 'Clusters');

Figura: Clústeres obtenidos a partir de la clasificación no supervisada

Asigna etiquetas a los clústeres

Tras una inspección visual, los clústeres obtenidos en los pasos anteriores coinciden estrechamente con los límites de la granja que se observan en la imagen de alta resolución. Según el conocimiento local, sabemos que los dos clústeres más grandes serían los de maíz y soja. Calculemos las áreas de cada clúster en nuestra imagen.

// Calculate Cluster Areas
// 1 Acre = 4046.86 Sq. Meters
var areaImage = ee.Image.pixelArea().divide(4046.86).addBands(clustered);

var areas = areaImage.reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'cluster',
    }),
    geometry: geometry,
    scale: 10,
    maxPixels: 1e10
    });

var clusterAreas = ee.List(areas.get('groups'));

// Process results to extract the areas and create a FeatureCollection

var clusterAreas = clusterAreas.map(function(item) {
  var areaDict = ee.Dictionary(item);
  var clusterNumber = areaDict.getNumber('cluster').format();
  var area = areaDict.getNumber('sum')
  return ee.Feature(null, {cluster: clusterNumber, area: area})
})

var clusterAreaFc = ee.FeatureCollection(clusterAreas);
print('Cluster Areas', clusterAreaFc);

Elegimos los 2 clústeres con el área más grande.

var selectedFc = clusterAreaFc.sort('area', false).limit(2);
print('Top 2 Clusters by Area', selectedFc);

Pero aún no sabemos qué clúster corresponde a qué cultivo. Si tuvieras algunas muestras de campo de maíz o soja, podrías superponerlas en los clústeres para determinar sus etiquetas respectivas. En ausencia de muestras de campo, podemos aprovechar las estadísticas agregadas de los cultivos. En muchas partes del mundo, se recopilan y publican periódicamente estadísticas agregadas sobre los cultivos. En el caso de EE.UU., el Servicio Nacional de Estadísticas Agrícolas (NASS) tiene estadísticas detalladas de los cultivos para cada condado y cada cultivo principal. En el año 2022, el condado de Cerro Gordo, Iowa, tuvo una superficie sembrada de maíz de 65,350 hectáreas y una superficie sembrada de soja de 44,720 hectáreas.

Con esta información, ahora sabemos que, entre los 2 clústeres principales, el que tiene el área más grande probablemente sea de maíz y el otro, de soja. Asignemos estas etiquetas y comparemos las áreas calculadas con las estadísticas publicadas.

var cornFeature = selectedFc.sort('area', false).first();
var soybeanFeature = selectedFc.sort('area').first();
var cornCluster = cornFeature.get('cluster');
var soybeanCluster = soybeanFeature.get('cluster');

print('Corn Area (Detected)', cornFeature.getNumber('area').round());
print('Corn Area (From Crop Statistics)', 163500);

print('Soybean Area (Detected)', soybeanFeature.getNumber('area').round());
print('Soybean Area (From Crop Statistics)', 110500);

Crea un mapa de cultivos

Ahora conocemos las etiquetas de cada clúster y podemos extraer los píxeles de cada tipo de cultivo y combinarlos para crear el mapa de cultivos final.

// Select the clusters to create the crop map
var corn = clustered.eq(ee.Number.parse(cornCluster));
var soybean = clustered.eq(ee.Number.parse(soybeanCluster));

var merged = corn.add(soybean.multiply(2));
var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(merged.clip(geometry), cropVis, 'Crop Map (Detected)');

Para ayudar a interpretar los resultados, también podemos usar elementos de la IU para crear y agregar una leyenda al mapa.

// Add a Legend
var legend = ui.Panel({
  layout: ui.Panel.Layout.Flow('horizontal'),
  style: {position: 'bottom-center', padding: '8px 15px'}});

var addItem = function(color, name) {
  var colorBox = ui.Label({
    style: {color: '#ffffff',
      backgroundColor: color,
      padding: '10px',
      margin: '0 4px 4px 0',
    }
  });
  var description = ui.Label({
    value: name,
    style: {
      margin: '0px 10px 0px 2px',
    }
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')}
)};

var title = ui.Label({
  value: 'Legend',
  style: {fontWeight: 'bold',
    fontSize: '16px',
    margin: '0px 10px 0px 4px'}});

legend.add(title);
legend.add(addItem('#ffd400', 'Corn'));
legend.add(addItem('#267300', 'Soybean'));
legend.add(addItem('#bdbdbd', 'Other Crops'));

Figura: Mapa de cultivos detectados con cultivos de maíz y soja

Valida los resultados

Pudimos obtener un mapa de tipos de cultivos con el conjunto de datos de Satellite Embedding sin etiquetas de campo, solo con las estadísticas agregadas y el conocimiento local de la región. Comparemos nuestros resultados con el mapa oficial de tipos de cultivos de las capas de datos de tierras de cultivo (CDL) del NASS del USDA.

var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var cropMap = cropLandcover.updateMask(croplandMask).rename('crops');

// Original data has unique values for each crop ranging from 0 to 254
var cropClasses = ee.List.sequence(0, 254);
// We remap all values as following
// Crop     | Source Value | Target Value
// Corn     | 1            | 1
// Soybean  | 5            | 2
// All other| 0-255        | 0
var targetClasses = ee.List.repeat(0, 255).set(1, 1).set(5, 2);
var cropMapReclass = cropMap.remap(cropClasses, targetClasses).rename('crops');

var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(cropMapReclass.clip(geometry), cropVis, 'Crop Landcover (CDL)');

Si bien existen discrepancias entre nuestros resultados y el mapa oficial, notarás que pudimos obtener resultados bastante buenos con un esfuerzo mínimo. Si aplicamos pasos de posprocesamiento a los resultados, podemos quitar parte del ruido y completar los vacíos en el resultado.

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