Modelo de detección de áreas verdes en ciudades — TUTORIAL Ciudad de Buenos Aires

En notas anteriores analizamos la importancia de la detección de áreas verdes para impulsar políticas públicas de calidad y cuidar el medioambiente. Además, revisamos los pasos principales para su mapeo a través del uso de herramientas de machine learning e imágenes satelitales: pre-procesamiento, entrenamiento, predicción y post-procesamiento.

En este posteo nos detendremos en detalle en cada una de esas etapas y te mostraremos paso a paso cómo detectar áreas verdes en ciudades y clasificarlas aplicando un modelo de clasificación basado en redes neuronales convolucionales.

Al final de este tutorial obtendremos un mapa de los espacios verdes de la Ciudad de Buenos Aires -nuestro caso de ejemplo-, clasificadas en dos clases: parques y vegetación (arbolado público, bulevares y terrenos baldíos).

Áreas verdes detectadas y clasificadas como parque (izquierda) y vegetación (derecha).

Descarga de datos

Comenzaremos descargando los datasets de interés. En este tutorial utilizaremos las herramientas de código libre Unetseg y Satproc, imágenes aéreas de la Ciudad de Buenos Aires y anotaciones de parques y vegetación de la ciudad.

• Descarga de imágenes aéreas de barrios de Buenos Aires (incluye Belgrano, Chacarita, Colegiales, Palermo, Recoleta, Villa Crespo).

import os
from glob import glob
barrio='Belgrano'
path_to_images= './workshop/imagenes/'
!mkdir -p $path_to_images #crea el directorio
bucket_imagen_caba=os.path.join('gs://dym-workshops-public/gcba/workshop/imagen_2021/2021_RGB_utm.tif') #imagen de CABA utilizada para entrenar
!gsutil -m cp -r $bucket_imagen_caba $path_to_images #descarga la imagen

• Descarga del archivo vectorial de anotaciones: este archivo contiene polígonos de las clases parques y vegetación.

path_to_gt='./workshop/data/gt/'
!mkdir -p $path_to_gt
bucket_gt='gs://dym-workshops-public/gcba/workshop/anotaciones/areas_verdes.gpkg'
!gsutil -m cp -r  $bucket_gt  $path_to_gt

• Descarga del vector de área de interés: el archivo vectorial de área de interés delimita las zonas donde están las anotaciones. Al definirlo, reducimos el tiempo de cómputo ya que no se recorre toda la imagen.

path_to_aoi='./workshop/data/shp/ '
!mkdir -p $path_to_aoi
bucket_aoi='gs://dym-workshops-public/gcba/workshop/data/shp/aoi_areas_verdes.geojson'
!gsutil -m cp -r  $bucket_aoi  $path_to_aoi

Pre-procesamiento: herramienta Satproc

El objetivo del pre-procesamiento es generar el dataset de entrenamiento con imágenes de buena calidad, cantidad y formato necesario. Los modelos de segmentación requieren un dataset conformado por imágenes y máscaras que delimitan el objeto de interés; estas son imágenes binarias con 1 donde está el objeto y 0 donde no.

Para ello vamos a ayudarnos con Satproc, nuestra herramienta que permite precisamente crear un dataset de imágenes y máscaras con el formato requerido por algunos modelos de machine learning de segmentación de imágenes.

¿Cómo funciona? Satproc recorre cada imagen con una ventana de búsqueda de un tamaño fijo (size). Si encuentra un polígono dentro, entonces genera una máscara con valores de pixel igual a 1 dónde está el polígono y 0 donde no. Luego se mueve un paso (step_size), y repite. Si hay más de una clase, entonces concatena las máscaras en una sola ventana formando una imagen de n-bandas, con n el número de clases.

Luego, toma la imagen y el archivo vectorial con las anotaciones y genera un directorio con dos carpetas: una es el dataset de imágenes y la otra el de máscaras (una máscara por cada imagen). Dado que queremos detectar dos clases de objetos, la máscara será una imagen binaria de 2 bandas, es decir, una para vegetación y otra para parques.

Ejemplo de una imagen (izquierda) y sus respectivas máscaras: máscara de vegetación (medio) y máscara de parque (derecha)

1.1 Generación del dataset de entrenamiento

Como mencionamos anteriormente, para la creación del dataset de entrenamiento necesitamos las imágenes satelitales con la combinación de bandas que deseamos utilizar y un archivo vectorial con la anotación de los polígonos de verdad de campo (.geojson .shp .gpkg etc), ambos con la misma georeferencia.

Para este estudio utilizamos EPSG:32721. Además, el archivo vectorial debe contener un campo “class” donde cada polígono esté identificado con su clase, como se muestra en la figura.

Class
1    P
2    V
3    V
4    V
5    P
6    V

>Instalar Satproc y algunas otras librerías que utilizaremos

!pip install pysatproc

> Para generar las imágenes del dataset de entrenamiento:

size = 3000
step_size = 1000
path_dataset_train=os.path.join("./workshop/dataset/data_train/" +str(size)+'_'+str(step_size)+'/')
!mkdir -p $path_dataset_train
vector_file = './workshop/data/gt/areas_verdes.gpkg'
aoi= "./workshop/data/shp//aoi_areas_verdes.geojson"
path_to_files = './workshop/imagenes/2021_RGB_utm.tif' #Imagen de CABA
output_folder = os.path.join("./workshop/dataset/data_train/" +str(size)+'_'+str(step_size)+'/')
!satproc_extract_chips \
$path_to_files \
-o  $output_folder \
--size $size \
--step-size $step_size \
--aoi $aoi \
--labels $vector_file \
--label-property 'class' \
--classes 'V' 'P'

Los argumentos son:

• o: es la ruta de destino > recomendamos que dicha ruta sea descriptiva, por ejemplo “data_train/3000_1000/ ” describe : Data_train → datos usados para entrenar; 3000_1000 → (las imágenes son cuadradas de 3000x3000 y el step_size es de 1000)
• size: tamaño de las imágenes resultantes (las imágenes son cuadradas)
• step-size: paso del proceso. Si step-size es igual que el size entonces no hay overlap en las imágenes.
• En ocasiones es útil para el entrenamiento generar los chips con un overlap de este modo tenemos más datos para entrenar. Pero en la predicción el valor debe ser igual al tamaño que la imagen.
• aoi: ruta al archivo vectorial donde están definidas los barrios. Al definir una región de interés sólo se procesan las imágenes que interceptan esos barrios.
• labels: es el archivo vectorial de anotaciones.
• label-property: nombre del campo donde se define cada categoría (solo se usa para el entrenamiento)
• classes: nombres de las clases (como aparecen en el .gpkg), separados por espacios

Este comando va a generar dos carpetas en la ruta de destino : “images” y “extent”. Los archivos de la primera van a ser de tipo Tiff de 3 bandas (rgb) y los de la segunda van a ser también de tipo Tiff, pero de N bandas donde N representa el número de clases. En este caso son 2 bandas y cada una de ellas es una máscara binaria.

2. Entrenamiento del modelo

Este modelo es una red neural profunda de segmentación de objetos de tipo CNN con una arquitectura U-Net. Originalmente, fue desarrollada para la segmentación de imágenes biomédicas pero se ha descubierto que tiene un buen desempeño en muchas otras áreas.

La arquitectura posee múltiples capas de down-sampling y up-sampling, las cuales están conectadas. De este modo se combina la información de capas que poseen activaciones de mayor precisión con los nodos de las capas anteriores que poseen mayor resolución.

Arquitectura U-net

• Input: dataset de imágenes y máscaras
• Output: por cada imagen el modelo devuelve una imagen donde el valor de cada píxel es la probabilidad de encontrar al objeto de interés. Si se predice sobre varias clases, entonces la imagen de salida será de n-bandas , una por cada clase.

Este tipo de análisis cuenta con la ventaja de que las imágenes pueden ser no solo “visuales” RGB. Se le pueden añadir otras bandas y también información georeferenciada como otras capas GIS para mejorar los resultados.

> Instalamos 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.

!pip install folium==0.2.1
!pip install unetseg
from unetseg.train import TrainConfig, train
from unetseg.evaluate import plot_data_generator
import os

> Configuramos el entrenamiento:

Obs: Es útil usar un nombre para el archivo de pesos (nombre del modelo) que dé información sobre los parámetros de entrenamiento. Por ejemplo: < dim de las imagenes > < num de bandas > < size >_< stepsize > < step_per_epoch > .h5 o similares

model_name='UNet_256x256_3N_2C_spe92_areas_verdes_modelo1.h5'

> Veamos cuántas imágenes y máscaras hay en el dataset de entrenamiento para poder calcular el valor del parámetro: steps_per_epoch.

from glob import glob
path_to_images=glob(os.path.join('../workshop/dataset/data_train/' +str(size)+'_'+str(step_size)+'/'+'images/*.tif'))
#path_to_mask=glob(os.path.join("../workshop/dataset/data_train/" +str(size)+'_'+str(step_size)+'/'+'extent/*.tif'))
n=(len(path_to_images),len(path_to_mask))
print(n)
config = TrainConfig(width=256,  #  tamaño de la imagen procesada por la UNet (debe ser multiplos de 16 , por ej 160, 320,etc; y no menor a 80)
height=256,
n_channels=3,  #  número de canales de la imagen, rgb -> 3
n_classes=2, # número de clases a clasificar
apply_image_augmentation=True, #  si es True , amplia el dataset generando imagenes nuevas a partir de pequeñas variaciones de las ya existentes (rotación,)
seed=42,
epochs=20, # Cantidad de veces que el dataset entero puede pasar por el proceso de entrenamiento
batch_size=16, #cantidad de datos que se procesan por vez, puede ser limitado por la memoria de gpu disponible (debe ser multiplo de 16)
steps_per_epoch=92, #  típicamente debe ser igual al numero de imágenes / el batch_size, si es mayor incrementara el número de imágenes generadas con image augmentation
early_stopping_patience=5, # a medida que entrena se guardan los resultados del entrenamiento despues de cada epoch, si el error no varió luego de ¿¿10 ?? iteraciones , se corta el proceso porque se entiende que el error ya disminuyó significativamente
validation_split=0.2, # se divide la muestra en training (80% de los datos) y validation (20% de los datos) para calcular el error durante el proceso de entrenamiento
test_split=0.1, # Cantidad de imágenes del dataset
images_path=os.path.join("../workshop/dataset/data_train/" +str(size)+'_'+str(step_size)+'/'),#ruta a las imágenes generadas con Satproc
                    model_path=os.path.join('../data/weights/', model_name),#  ruta del modelo entrenado
model_architecture='unet',
evaluate=True,
class_weights= [0.2, 0.8])

plot_data_generator permite visualizar algunas de las imágenes y máscaras del dataset de entrenamiento.

plot_data_generator(num_samples=3, fig_size=(10, 10), train_config=config,img_ch = 2)

> Corremos el entrenamiento

res_config = train(config)

Métrica del modelo

Una de las métricas más utilizadas en computer vision es Intersection over Unión (IoU).

Ésta calcula el área de intersección entre la verdad de campo y la predicción sobre el área común a ambos. Tiende a 1 a medida que la predicción mejora.

Cáculo de la Intersection over Unión (IoU).

import matplotlib.pyplot as plt
plt.figure(figsize=(16,4))
plt.subplot(121)
plt.plot(res_config.history['loss'])
plt.plot(res_config.history['val_loss'])
plt.title('Loss')
plt.ylabel('Loss')
plt.xlabel('Epoch')
plt.legend(['Train', 'Val'], loc='upper left')
plt.subplot(122)
plt.plot(res_config.history['mean_iou'])
plt.plot(res_config.history['val_mean_iou'])
plt.title('mean_iou')
plt.ylabel('val_mean_iou')
plt.xlabel('Epoch')
plt.legend(['Train', 'Val'], loc='upper left')
plt.show()

3) Predicción

Ahora crearemos el dataset de predicción. Con él, el modelo que entrenamos previamente generará las predicciones sobre estas nuevas imágenes devolviendo la probabilidad de encontrar los objetos de interés. Para facilitar este proceso utilizaremos Unetseg.

!pip install folium==0.2.1
!pip uninstall -y h5py
!pip install h5py==2.10.0
!pip install unetseg
!pip install pysatproc
from unetseg.predict import PredictConfig, predict
from unetseg.evaluate import plot_data_results
import os

> Descargamos la imagen del barrio

barrio='Belgrano'
path_to_images=os.path.join('./workshop/images_predict/'+barrio+'/')
!mkdir -p $path_to_images
bucket_imagen_barrio=os.path.join('gs://dym-workshops-public/gcba/workshop/'+barrio+'/imagen/')
!gsutil -m cp -r $bucket_imagen_barrio $path_to_images

> Descargamos el archivo vectorial de área de interés

bucket_aoi_barrio=os.path.join('gs://dym-workshops-public/gcba/workshop/'+barrio+'/shp/'+barrio+'_1.geojson')
aoi_predict='./workshop/aoi_predict/'
!mkdir $aoi_predict
!gsutil -m cp -r $bucket_aoi_barrio $aoi_predict

> Descargamos el modelo ya entrenado

path_to_model='./workshop/data/'
!mkdir -p path_to_model
bucket_modelo='gs://dym-workshops-public/gcba/workshop/modelo'
!gsutil -m cp -r $bucket_modelo $path_to_model

Ahora sí, generamos el dataset de predicción. De modo similar, utilizamos satproc_extract_chips para generar el dataset de predicción. En este caso el archivo vectorial delimita las zonas de interés para la predicción y el size y step_size son iguales, ya que no es útil un overlap de las imágenes.

size = 3000
step_size=size
path_to_files = os.path.join('./workshop/images_predict/'+barrio+'/imagen/*.tif') #imagen del barrio
output_folder = os.path.join('./workshop/'+barrio+'/dataset/data_predict/'+str(size)+'_'+str(step_size)+'/')
vector_file_aoi = os.path.join('./workshop/aoi_predict/'+barrio+'_1.geojson')
!satproc_extract_chips \
$path_to_files \
-o  $output_folder \
--size $size \
--step-size $step_size \
--aoi $vector_file_aoi

> Consultamos el nombre del modelo

!ls ./workshop/data/

Definimos la configuración para la predicción. Debemos pasar la ruta a las imágenes del dataset de predicción y al modelo.

predict_config = PredictConfig(
images_path=os.path.join('./workshop/'+barrio+'/dataset/data_predict/'+str(size)+'_'+str(step_size)+'/'), # ruta a las imagenes sobre las cuales queremos predecir
results_path=os.path.join('./workshop/'+barrio+'/dataset/data_results/'+str(size)+'_'+str(step_size)+'/'), # ruta de destino para nuestra predicción
batch_size=16,
model_path=os.path.join('./workshop/data/','UNet_256x256_3N_2C_spe92_areas_verdes_modelo1.h5'),  #  ruta al modelo (.h5)
height=256,
width=256,
n_channels=3,
n_classes=2,
class_weights= [0.2,0.8])
predict(predict_config)  # Ejecuta la predicción
plot_data_results(num_samples=2, fig_size=(5, 5), predict_config=predict_config, img_ch =2,n_bands =3)

4) Post-procesamiento de los resultados

Finalmente llegamos al post-procesamiento, que tiene como objetivo mejorar la calidad de los resultados. En esta etapa se aplica un umbral para conservar aquellas predicciones de mayor probabilidad. Asimismo, una vez poligonizado el resultado final, se puede aplicar un filtro de mínima área para disminuir la cantidad de predicciones “ruidosas”.

La siguiente función aplica una rutina de poligonización sobre los resultados de la predicción del modelo y genera un archivo vectorial en formato GeoPackage. (GPKG). La rutina aplica gdal_polygonize.py a cada chip resultado generando un GPKG para cada chip, y luego une todos estos archivos en uno solo, de manera eficiente.

Antes de unirlos también aplica un umbral sobre los valores de los rásteres, que en este caso representan la probabilidad (valores entre 0 y 1). filter_by_max_prob es una función de Satproc que filtra los resultados donde ningún pixel supera el umbral definido.

> Ahora, definimos un umbral de probabilidad para cada clase:

threshold_A=50
threshold_B=50

El siguiente código aplica un umbral a la probabilidad para cada clase y asigna:

• Valor 100 a los píxeles de la clase “Parque”
• Valor 200 a los píxeles de la clase “Vegetación”
• Valor 1 al background

Esto se hace con el fin de transformar la imagen de probabilidades de 2 bandas (binarias), en una imagen de una banda con dichos valores de píxeles.

from tqdm import tqdm
import subprocess
import os
DEBUG=True
results_path = os.path.join('./workshop/'+barrio+'/dataset/data_results/'+str(size)+'_'+str(step_size)+'/') #ruta de los resultados de la predicción
#if DEBUG: print (results_path)
dest = os.path.join('./workshop/'+barrio+'/dataset/data_results/temp_postprocesamiento',str(size)+"_"+str(step_size)) # ruta de destino para el post-procesamiento
os.makedirs(dest, exist_ok=True)
#if DEBUG: print (dest)
files = next(os.walk(results_path))[2] #recorre los archivos
new_dir = os.path.join(dest, '1band_2classes','thr_'+str(threshold_A)+'_'+str(threshold_B)) #creamos un directorio para la imagen de 1 banda
os.makedirs(new_dir, exist_ok=True) #si no existe el directorio, lo crea
#if DEBUG: print (new_dir)
for file in files:
if os.path.splitext(file)[-1]=='xml': continue
result_dst = os.path.join(results_path,file)
output_tif_path = os.path.join(dest, '1band_2classes/','thr_'+str(threshold_A)+'_'+str(threshold_B),file) #destino de la imagen de una banda
exp = "((B<"+str(threshold_B)+")*(A >"+str(threshold_A)+")*199 ) + ((A<"+str(threshold_A)+")*(B >"+str(threshold_B)+")*99 ) +1" #cálculo para aplicar thresholds y asignar valores a los píxeles segun la clase
cmd_calc = f'gdal_calc.py -A {result_dst} --A_band=1 -B {result_dst} --B_band=2  --outfile {output_tif_path} --calc="{exp}" --NoDataValue=0'
subprocess.run(cmd_calc, shell=True)
print("--",output_tif_path)

> Transformamos la imagen de probabilidades de 2 bandas, en una imagen de una banda:

output_vrt = os.path.join(dest, ‘1band_2classes’,’thr_’+str(threshold_A)+’_’+str(threshold_B),’tiles.vrt’)
output_tif = os.path.join(dest, ‘1band_2classes’,’thr_’+str(threshold_A)+’_’+str(threshold_B),’*.tif’)
cmd_buildvrt = f’gdalbuildvrt {output_vrt} {output_tif}’
subprocess.run(cmd_buildvrt, shell=True)

> Realizamos el poligonizado:

output_shp_path = os.path.join('./workshp/'+barrio+'/dataset/data_results/', 'poligonize_classAB_thr_'+str(threshold_A)+'_'+str(threshold_B)+'_'+str(size)+'_'+str(step_size)+'.gpkg')
cmd_polygonize = f'gdal_polygonize.py {output_vrt} {output_shp_path}'
subprocess.run(cmd_polygonize, shell=True)

De manera opcional, podemos aplicar un filtro de mínima área primero debemos pasar los resultados a metros (utm):

#src_file = 'file_4326.gpkg'
#dst_file = 'file_utm_32721.gpkg'
#!ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:32721 -f 'GPKG' $dst_file $src_file

> Finalmente, filtramos aquellos polígonos que no cumplen con el valor de área mínima definido:

min_area = 10
input_path = dst_file
output_path = "file_utm_32721.gpkg"
!ogr2ogr \
-t_srs EPSG:32721 \
-f "GPKG" \
-sql "SELECT * FROM tabla m WHERE (ST_Area(geom) > $min_area) " \
-dialect SQLITE \
-nln results \
$output_path \
$input_path

Si seguiste todos estos pasos, habrás obtenido imágenes como esta:

Detección y clasificación de áreas verdes en vegetación (rojo) y parques (azul)

Te dejamos el link al Google Colab para que puedas correr el código desde allí y acceder a más información sobre la detección de áreas verdes en ciudades.

Si tienes dudas o quieres saber más sobre nuestras soluciones te invitamos a seguirnos por nuestras redes sociales o escribirnos.