Cómo crear un modelo de cobertura del suelo en Sudamérica en 4 pasos

Recientemente, Radiant Earth Foundation publicó un conjunto de datos de cobertura de suelos para Sudamérica, continuando la labor que vienen realizando en otras partes del mundo y áreas de interés.

Siguiendo el hilo de publicaciones anteriores, en este posteo mostraremos cómo se puede entrenar un modelo de segmentación a partir de este conjunto de datos y en solo 4 pasos. Puntualmente, vamos a detallar cómo entrenar un modelo para la clasificación del uso del suelo asociado a cultivos, a partir del conjunto de datos mencionado.

El dataset: LandCoverNet South America

El dataset publicado comprende un conjunto de datos de etiquetas e imágenes satelitales de las misiones Sentinel-1, Sentinel-2 y Landsat 8 para clasificación de usos del suelo del área sudamericana (si quieres saber más sobre fuentes de imágenes satelitales, ingresa aquí).

Cada pixel es identificado dentro de las siete posibles clases de tipo de superficie: de agua, suelo desnudo natural, suelo desnudo artificial, vegetación leñosa, suelo cultivado, suelo semi-cultivados y suelo de hielo o nieve permanente.

Esta primera versión de LandCoverNet South America ofrece un total de 1200 chips de imágenes de 256 x 256 píxeles, y abarca 40 mosaicos. En cada chip encontraremos observaciones temporales con una etiqueta de clasificación anual en formato raster (GeoTIFF files) de los siguientes satélites:

• Sentinel-1: detección GRD (ground range distance) con calibración radiométrica y ortorrectificación a una resolución espacial de 10 m.
• Sentinel-2: producto de reflectancia superficial (L2A) a una resolución espacial de 10 m
• Landsat-8: producto de reflectancia superficial de Collection 2 Level-2

Entre las ventajas del uso de LandCoverNet, encontramos que su licencia CC-BY-4.0 permite la utilización y modificación de su contenido sin grandes restricciones. Para conocer más sobre este desarrollo de Radiant Earth, puedes ingresar aquí*.

Entrenar un modelo para la clasificación del uso del suelo asociado a cultivos

Antes de comenzar, se sugiere enfáticamente ejecutar este código en un entorno que permita acceder a una GPU. Al final de esta guía se deja el link para acceder al notebook en formato Google Colaboratory.

Pasos:

1- Descarga de los conjuntos de datos

2- Preprocesamiento imágenes sentinel-2

3- Preprocesamiento anotaciones

4- Entrenamiento del modelo

1- Descarga del conjunto de datos

Para avanzar con este paso, sugerimos ver la guía publicada por Radiant Earth. Incluye la registración como usuario y la generación de un API. Los datos.

import os
os.environ['MLHUB_API_KEY'] = <API KEY>

En esta guía vamos a trabajar con imágenes Sentinel-2 y las anotaciones. En particular, descargaremos las bandas asociadas al Azul, Verde y Rojo (2, 3 y 4). Limitaremos a cien casos (en lugar de los 1.200 disponibles) para acelerar el proceso. Para obtener datos relacionados a la licencia, cita y descripción, ejecutamos:

collection_id = "ref_landcovernet_sa_v1_source_sentinel_2"
collection = client.get_collection(collection_id)
print(f'Description: {collection["description"]}')
print(f'License: {collection["license"]}')
print(f'DOI: {collection["sci:doi"]}')
print(f'Citation: {collection["sci:citation"]}')
Description: LandCoverNet South America Sentinel 2 Source Imagery
License: CC-BY-4.0
DOI: doi.org/10.34911/rdnt.6a27yv
Citation: Radiant Earth Foundation (2022) "LandCoverNet South America: A Geographically Diverse Land Cover Classification Training Dataset", Version 1.0, Radiant MLHub. https://doi.org/10.34911/rdnt.6a27yv

Para descargarlos, vamos a usar un conjunto de funciones que vienen como ejemplo en la documentación. Si bien no las agregamos a este post para que no sea tan largo, en el notebook final las incluimos y también las pueden encontrar en la documentación oficial. Ejecutamos:

items = get_items(collection_id, max_items=100)
for item in items:
    download_labels_and_source(item, assets=['B02', 'B03', 'B04'])

Los chips que se descargaron se almacenan por defecto en el directorio data.

Para el caso de las anotaciones, descargamos directamente el archivo fuente (comprimido).

collection_id = 'ref_landcovernet_sa_v1_labels'
collection = client.get_collection(collection_id)
print(f'Description: {collection["description"]}')
print(f'License: {collection["license"]}')
print(f'DOI: {collection["sci:doi"]}')
print(f'Citation: {collection["sci:citation"]}')
Description: LandCoverNet South America Labels
License: CC-BY-4.0
DOI: doi.org/10.34911/rdnt.6a27yv
Citation: Radiant Earth Foundation (2022) "LandCoverNet South America: A Geographically Diverse Land Cover Classification Training Dataset", Version 1.0, Radiant MLHub. https://doi.org/10.34911/rdnt.6a27yv

Descargamos el archivo comprimido en el directorio labels. Luego lo descomprimimos allí.

client.download_archive(collection_id, output_dir='./labels')
!tar -xvzf ./labels/ref_landcovernet_sa_v1_labels.tar.gz -C ./labels/

El objetivo de los siguientes pasos de preprocesamiento viene dado por acomodar el formato que posee el conjunto de datos descargado a la estructura de directorios que requiere el paquete unetseg. Refrescando este punto, se requiere que los chips con las imágenes se encuentren en un directorio llamado images, mientras que las anotaciones tienen que estar en el directorio extent. El nombre de cada chip debe ser el mismo en cada directorio.

Puede sonar raro para aquellos que vengan siguiendo nuestros posteos el hecho de que no usemos el paquete satproc aquí. Esto se debe a que los datos ya vienen en formato chip. Sin embargo, si se requiriera cambiar el tamaño del chip, satproc podría ser una buena opción para realizar ese paso.

2- Preprocesamiento de imágenes Sentinel-2

La estructura que tiene cada chip, a partir de haber descargado las bandas 2, 3 y 4, es:

!ls data/ref_landcovernet_sa_v1_source_sentinel_2_21MXR_18_20181231/labels
21MXR_18_20181231_B02_10m.tif  21MXR_18_20181231_B04_10m.tif
21MXR_18_20181231_B03_10m.tif

A partir de la mencionada estructura de directorios necesaria para entrenar con unetseg, necesitaremos entonces renombrar cada archivo raster con el nombre del chip correspondiente, para luego moverlos al directorio images. Previamente, vamos a unir los tres archivos asociados a cada banda en uno solo, utilizando un raster virtual (que luego borramos).

chips = glob.glob('./data/ref_landcovernet_sa_v1_source_sentinel_2*')
os.makedirs('./train/images', exist_ok=True)
for chip_dir in chips:
    chip_name_vrt = chip_dir.replace('./data/ref_landcovernet_sa_v1_source_sentinel_2_','') + '.vrt'
    chip_name_tif = ('_').join(chip_dir.replace('./data/ref_landcovernet_sa_v1_source_sentinel_2_','').split('_')[:2])+ '.tif'
    print(chip_name_tif)
    outfile_vrt = os.path.join('train/images', chip_name_vrt)
    outfile_tif = os.path.join('train/images', chip_name_tif)
    inputs_path = os.path.join(chip_dir, 'labels', '*.tif')
!gdalbuildvrt -separate $outfile_vrt $inputs_path
!gdal_translate $outfile_vrt $outfile_tif
!rm $outfile_vrt

La estructura resultantes es:

!ls images
21LUE_25.tif  21MXR_15.tif  23LKC_00.tif  23LKC_20.tif 23LMJ_10.tif
21LUE_26.tif  21MXR_16.tif  23LKC_01.tif  23LKC_21.tif 23LMJ_11.tif
21LUE_27.tif  21MXR_17.tif  23LKC_02.tif  23LKC_22.tif 23LMJ_12.tif
21LUE_28.tif  21MXR_18.tif  23LKC_03.tif  23LKC_23.tif 23LMJ_13.tif
21LUE_29.tif  21MXR_19.tif  23LKC_04.tif  23LKC_24.tif 23LMJ_14.tif
21MXR_00.tif  21MXR_20.tif  23LKC_05.tif  23LKC_25.tif 23LMJ_15.tif
21MXR_01.tif  21MXR_21.tif  23LKC_06.tif  23LKC_26.tif 23LMJ_16.tif
21MXR_02.tif  21MXR_22.tif  23LKC_07.tif  23LKC_27.tif 23LMJ_17.tif
21MXR_03.tif  21MXR_23.tif  23LKC_08.tif  23LKC_28.tif 23LMJ_18.tif
21MXR_04.tif  21MXR_24.tif  23LKC_09.tif  23LKC_29.tif 23LMJ_19.tif
21MXR_05.tif  21MXR_25.tif  23LKC_10.tif  23LMJ_00.tif 23LMJ_20.tif
21MXR_06.tif  21MXR_26.tif  23LKC_11.tif  23LMJ_01.tif 23LMJ_21.tif
21MXR_07.tif  21MXR_27.tif  23LKC_12.tif  23LMJ_02.tif 23LMJ_22.tif
21MXR_08.tif  21MXR_28.tif  23LKC_13.tif  23LMJ_03.tif 23LMJ_23.tif
21MXR_09.tif  21MXR_29.tif  23LKC_14.tif  23LMJ_04.tif 23LMJ_24.tif
21MXR_10.tif  22KEU_05.tif  23LKC_15.tif  23LMJ_05.tif 23LMJ_25.tif
21MXR_11.tif  22KEU_07.tif  23LKC_16.tif  23LMJ_06.tif 23LMJ_26.tif
21MXR_12.tif  22KEU_13.tif  23LKC_17.tif  23LMJ_07.tif 23LMJ_27.tif
21MXR_13.tif  22KEU_21.tif  23LKC_18.tif  23LMJ_08.tif 23LMJ_28.tif
21MXR_14.tif  22KEU_22.tif  23LKC_19.tif  23LMJ_09.tif 23LMJ_29.tif

3- Preprocesamiento de anotaciones

En este caso, la estructura que tiene cada chip es la siguiente:

ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/labels.tif
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/source_dates.csv
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/stac.json

El archivo raster con las anotaciones viene en el labels.tif. Cabe destacar que los chips corresponden a imágenes de diferentes fechas. En este posteo ignoraremos ese hecho, considerando, para simplificar el análisis, a todos los chips como si fueran de la misma fecha. Cabe destacar que el valor de píxel asociado a la clase de Áreas Cultivada es 6 (lo usaremos en el comando gdal_calc.py).

os.makedirs('extent', exist_ok=True)
chips = glob.glob('data/ref_landcovernet_sa_v1_labels/*/*.tif')
for chip_path in chips:
    chip_dir = os.path.dirname(chip_path)
    chip_name = chip_dir.replace('data/ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_','')
    chip_name_tif = os.path.join('./train/extent',chip_name + '.tif')
!gdal_calc.py -A $chip_path --outfile=$chip_name_tif --calc="A==6" --NoDataValue=0

La estructura resultantes es:

ref_landcovernet_sa_v1_labels/
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/labels.tif
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/source_dates.csv
ref_landcovernet_sa_v1_labels/ref_landcovernet_sa_v1_labels_21MXR_18/stac.json

Ahora que tenemos los archivos estructurados, tanto en su nombre como en los directorios que ocupan, procedemos a configurar el entrenamiento del modelo.

4- Entrenamiento del modelo

En primer lugar, instalamos el paquete de python unetseg. Este paquete sirve para simplificar el entrenamiento y predicción de modelos de deep learning, basados en la arquitectura U-Net, para problemas de segmentación semántica sobre imágenes geoespaciales.

!pip3 install unetseg
from unetseg.train import TrainConfig, train
from unetseg.evaluate import plot_data_generator

Luego, instanciamos el entrenamiento, configurando los parámetros disponibles. En la documentación se listan aquellos que actualmente se ofrecen para entrenar los modelos.

train_config = TrainConfig(
    width=256,
    height=256,
    n_channels=3,
    n_classes=1,
    apply_image_augmentation=True,
    seed=42,
    epochs=15,
    batch_size=4,
    steps_per_epoch=1,
    early_stopping_patience=3,
    validation_split=0.2,
    test_split=0.1,
    model_architecture="unet",
    images_path=os.path.join('./model'),
    model_path=os.path.join('./model/weights/', 'UNet_test.h5'),
    evaluate=True,
)

Como se ve allí, el tamaño de los chips es de 256x256 pixeles, tal como vienen los chips provistos en el dataset. La cantidad de clases es 1, tal como se planteó en la sección anterior. Los pesos del modelo, que usaremos para predecir en el siguiente paso, se guardan en el directorio model. El resto de los parámetros se seleccionaron aleatoriamente, sugiriendo al usuario experimentar con ellos para obtener la mejor combinación.

Veamos un ejemplo de cómo lucen los chips y las máscaras generadas anteriormente:

plot_data_generator(train_config=train_config, num_samples=1, fig_size=(10,5))

Entrenamos el modelo:

res_config = train(train_config)

Una vez entrenado el modelo, este se puede utilizar para predecir en nuevas imágenes. En un siguiente posteo ahondaremos sobre esto. Para aquellos ansiosos que quieran seguir con este paso, les sugerimos que vean el este tutorial.

Notebook: https://colab.research.google.com/drive/1QweySv8JCw1hj89nxt0SnVaPG5HuHPGY?usp=sharing

Si tienes dudas o interés en aplicar el código en tu proyecto o negocio, no dudes en escribirnos a contact@dymaxionlabs.com.

*Radiant Earth Foundation (2022) “LandCoverNet South America: A Geographically Diverse Land Cover Classification Training Dataset”, Version 1.0, Radiant MLHub. [Date Accessed] https://doi.org/10.34911/rdnt.6a27yv