Ejercicio: OBIA en GRASS GIS
En este ejercicio vamos a ejemplificar el flujo de trabajo para realizar una clasificación supervisada basada en objetos con datos SPOT.
Antes de empezar y para ganar tiempo, conectamos nuestro drive e instalamos GRASS en Google Colab.
Chequeamos el path de instalación de GRASS
e instalamos otras librerías de python que nos serán útiles.
La función que usaremos en esta sesión para realizar clasificación supervisada de los objetos llama internamente a ciertas librerías de R. Para poder tenerlas de antemano en nuestro entorno, primero instalamos y cargamos en la notebook el paquete de python que nos permite hacer interfaz con R.
Chequeamos nuestra sesión de R.
Instalamos todos los paquetes necesarios para este ejercicio. This might take a while…
Datos para esta sesión
- SPOT 6
- Canales VIS - NIR (6 m)
- Canal PAN (1.5 m)
- Datos ya corregidos y fusionados; resolución 1.5 m
Manos a la obra
Iniciamos GRASS GIS
Iniciamos GRASS GIS en el proyecto posgar2007_4_cba y mapset PERMANENT.
- Crear un mapset llamado
obia_spote importar allí las bandas VIS-NIR y PAN de la imagen SPOT con reproyección al vuelo y resolución espacial de 1.5 m. - Alinear la región computacional a la extensión y resolución de alguna de las bandas importadas.
- Hacer una ecualización de histograma para mejorar el contraste de visualización.
- Mostrar la combinación RGB color natural (1: azul, 2: verde, 3: rojo, 4: NIR).
Importar datos y visualizar
Creamos el mapset obia_spot.
Una vez creado el mapset, procedemos con la importación de datos. Empezamos por las bandas multi-espectrales.
Importamos también la banda pancromática.
Chequeamos la lista de mapas raster en el mapset.
Establecemos la región computacional a una de las bandas importadas, de modo que límites y resolución de la región queden alineados a los datos importados.
Establecemos grey como paleta de colores para las bandas RGB, para poder realizar la composición color natural posteriormente.
Primero, hacemos la ecualización de colores para mejorar la visualización.
Visualizamos la composición RGB 321 color natural.
Hay valores nulos?
En la cadena de procesamiento que vamos a realizar, especialmente en los pasos de segmentación, no se admiten valores nulos en los mapas de entrada. Así que, como primer paso antes de proceder, chequeamos la estadística univariada para saber si hay valores nuelos en nuestros mapas.
Si hubiera valores nulos, se deben rellenar antes de comenzar! Qué herramientas podríamos utilizar?
Índices espectrales y texturas GLCM
Como ya hicimos en los ejercicios anteriores, estimamos índices espectrales de vegetación y agua.
# install i.wi
gs.run_command("g.extension",
extension="i.wi")
# estimate water index
gs.run_command("i.wi",
green="SPOT_20180621_PANSHARP.2",
nir="SPOT_20180621_PANSHARP.4",
output="SPOT_20180621_NDWI",
winame="ndwi_mf")
# set ndwi color palette
gs.run_command("r.colors",
map="SPOT_20180621_NDWI",
color="ndwi")Por otra parte, vamos a estimar medidas de textura con r.texture. En este caso, elegimos Inverse Difference Moment (IDM) y Angular Second Moment (ASM), pero hay muchas otras opciones. Este módulo puede ser muy lento si la región computacional es muy grande y si el tamaño de la ventana es grande también. Para agilizar los cálculos en estos casos, se puede usar la extensión r.texture.tiled.
Visualizamos las bandas creadas. Esto puede no funcionar dentro de Colab, en tal caso probar de a un mapa o usar la clase Map.
idx_map = gj.InteractiveMap(width = 400, use_region=False, tiles="OpenStreetMap")
idx_map.add_raster("SPOT_20180621_NDVI", opacity=0.7)
idx_map.add_raster("SPOT_20180621_NDWI", opacity=0.7)
idx_map.add_raster("SPOT_20180621_IDM", opacity=0.7)
idx_map.add_raster("SPOT_20180621_ASM", opacity=0.7)
idx_map.add_layer_control(position = "bottomright")
idx_map.show()Sobre qué banda calculamos las texturas? Si no contamos con una banda pancromática, podemos crearla promediando las bandas visibles
Segmentación
Búsqueda de umbrales de sub y sobre-segmentación
Como vimos en el ejercicio anterior, usualmente, los módulos i.* toman un grupo de mapas como entrada. Antes de comenzar con la segmentación, entonces, creamos un grupo con las bandas multiespectrales únicamente.
Ahora vamos a aprovechar una de las grandes ventajas de la región computacional, i.e., definir una región más pequeña para realizar unas pruebas :)
Ahora sí, ejecutamos un par de segmentaciones para determinar niveles de sub- y sobre-segmentación. Empezamos con un umbral pequeño.
Probamos con un umbral más grande.
Hagamos zoom sobre los resultados para cada umbral.
Sobre-segmentado
Sub-segmentado
Se animan a probar con otros valores y en otras regiones?
Búsqueda automática de umbrales por optimización
Teniendo valores de umbrales de sub- y sobre-segmentación, podemos proceder a realizar la segmentación de manera autómatica, buscando el mejor umbral en una especie de grilla definida por umbral mínimo, umbral máximo y paso. Para ello, utilizamos la extensión: i.segment.uspo.
Este procedimiento es computacionalmente intensivo para un área grande y muchas combinaciones de parámetros, pero siempre podemos:
- Limitar el tamaño de la región computacional o generar tiles con i.cutlines y paralelizar la USPO.
- Limitar el rango de los parámetros, o usar un paso relativamente grande.
- Crear superpixels con i.superpixels.slic para usarlos como semillas.
Generación de semillas
La extensión i.superpixels.slic
- También puede utilizarse para la segmentación per se, como vimos en el ejercicio anterior.
- Es muy rápida para reagrupar pequeñas cantidades de píxeles similares.
- Se puede usar para reducir el número de píxeles en un factor de 4-5 y acelerar i.segment.uspo, que ya no empieza por pixeles individuales, sino por pequeños grupos de pixeles similares.
- Se usa baja compactación para mantener la separación espectral (Ver el manual: A larger compactness value will cause spatially more compact, but spectrally more heterogeneous superpixels).
USPO con superpixels como semillas
Vamos a ejecutar entonces i.superpixels.slic con bajo valor de compactación para generar el semillero para i.segment.uspo.
Acá podemos ver un recorte de la composición RGB y el resultado de la ejecución de i.superpixels.slic.
Instalamos las extensiones y ejecutamos la segmentación con optimización.
Convertimos el mejor resultado, i.e., **rank1* a vector para visualizar el resultado.
Zoom al resultado de ejecutar la segmentación con USPO
Cuántos segmentos obtuvieron?
Dar una mirada a v.info para responder.
Estadísticas de segmentos: i.segment.stats
Instalamos la extensión que nos permite extraer estadísticas para los segmentos obtenidos en el paso anterior.
# extract stats for all segments
# Note: *vectormap* output does not work for this dataset,
# see workaround with csv output in further steps
gs.run_command("i.segment.stats",
map="segs_obia_subset_rank1",
csvfile="segs_stats.csv",
rasters="SPOT_20180621_ASM,SPOT_20180621_IDM,SPOT_20180621_NDVI,SPOT_20180621_NDWI,SPOT_20180621_PAN",
raster_statistics="mean,stddev",
area_measures="area,perimeter,compact_circle,compact_square",
processes=4)Visualizamos el csv resultante con la estadística por segmentos.
Qué otras estadísticas se podrían obtener? Qué otro(s) módulo(s) podría(n) sustituir a i.segment.stats?
Datos de entrenamiento
Tenemos nuestros segmentos caracterizados, pero nos falta un dato fundamental para poder realizar una clasificación supervisada… la verdad de terreno, i.e., a qué clase pertenece una muestra de todos los segmentos del área de estudio. Como etiquetar segmentos o puntos es una tarea tediosa, se provee dentro del mapset PERMANENT un conjunto de puntos con etiqueta para realizar el entrenamiento. Veamos cuántos puntos tenemos.
Como el acceso a tablas de otros mapsets no está permitido, nos copiamos el vector con los puntos etiquetados al mapset obia_spot.
Cuántos puntos de cada clase tenemos?
Seleccionar segmentos sobre los cuales tenemos puntos de entrenamiento
Cuántos segmentos contienen puntos de entrenamiento?
Veamos un zoom a la selección de segmentos con puntos de entrenamiento.
Antes de asignar el atributo de los puntos a los segmentos, necesitamos agregar una tabla y una columna al vector con los segmentos.
Ahora sí, usamos el módulo v.distance para transferir la etiqueta (clase de cobertura) desde los puntos a los segmentos
Cuántos segmentos de cada clase tenemos?
Como, por alguna razón, para estos datos, no pudimos obtener un vector a partir de i.segment.stats, necesitamos hacer un “workaround” para obtener, en cambio, un csv con los segmentos de entrenamiento, sus estadísticas y sus etiquetas. Para eso, haremos lo siguiente:
- convertir a raster el vector
train_segmentscon las etiquetas, - obtener las estadísticas para los segmentos de entrenamiento con i.segment.stats,
- importar a GRASS la tabla csv resultante,
- unir el vector de segmentos de entrenamiento que tiene las etiquetas con la tabla con las estadísticas,
- eliminar la columna cat_ resultante de la unión, ya que no es un atributo para la clasificación, y
- exportar la nueva tabla de atributos como csv.
Empecemos el workaround…
# extract stats for training segments
gs.run_command("i.segment.stats",
map="train_segments",
csvfile="train_segs_stats.csv",
separator="comma",
rasters="SPOT_20180621_ASM,SPOT_20180621_IDM,SPOT_20180621_NDVI,SPOT_20180621_NDWI,SPOT_20180621_PAN",
raster_statistics="mean,stddev",
area_measures="area,perimeter,compact_circle,compact_square",
processes=4)Si de i.segment.stats ustedes pueden obtener un vector, la tarea es mucho más sencilla:
# select segments that are below labeled points
v.select \
ainput=segs_stats \ # vector obtenido con i.segment.stats
binput=labeled_points \
output=train_segments \
operator=overlap
# add column to train segments
v.db.addcolumn train_segments column="class int"
# assign label from points to segments
v.distance from=train_segments \
to=labeled_points \
upload=to_attr \
column=class \
to_column=train_class- Ejecutar una clasificación no supervisada con 10 clases
- Extraer una x cantidad de puntos por clase (r.sample.category)
- Etiquetar los puntos manualmente
- Usar puntos para transferir las etiquetas a los segmentos como ya vimos
# Unsupervised classification
i.group group=spot_all input=SPOT_20180621_ASM,SPOT_20180621_IDM,SPOT_20180621_NDVI,SPOT_20180621_NDWI,SPOT_20180621_PAN,SPOT_20180621_PANSHARP.1,SPOT_20180621_PANSHARP.2,SPOT_20180621_PANSHARP.3,SPOT_20180621_PANSHARP.4
i.cluster group=spot_all signaturefile=sig classes=10
i.maxlik group=spot_all signaturefile=sig output=uns_clas
# install extension
g.extension r.sample.category
# get n points per class
r.sample.category input=uns_clas output=uns_clas_points npoints=150
# then, manually label pointsClasificación con Machine learning
Ahora sí, finalmente llegamos a la clasificación de los segmentos. Vamos a instalar la extensión v.class.mlR. Esta extensión usa paquetes de R para realizar la clasificación, por lo tanto necesitamos tener R y ciertos paquetes previamente instalados.
# run classification
gs.run_command("v.class.mlR",
segments_file="segs_stats.csv", # stats of all segments
training_file="train_segs_stats_class.csv", # stats of training segs
train_class_column="class",
classified_map="classification",
raster_segments_map="segs_obia_subset_rank1",
classifier="rf",
folds=5,
partitions=10,
tunelength=10,
weighting_modes="smv",
weighting_metric="accuracy",
output_model_file="model",
variable_importance_file="var_imp.txt",
accuracy_file="accuracy.csv",
model_details="classifier_runs.txt",
r_script_file="Rscript_mlR.R",
processes=2, # use more if available
flags="n")Crear paleta de colores para las clases de cobertura urbana y guardarla como obia_urban en el directorio de trabajo. Ver el manual para obtener ejemplos sobre la creación de reglas. En la imagen de abajo, se sugieren algunos colores.
Aplicamos la paleta de colores que creamos.
Visualizamos el resultado de la clasificación supervisada con Machine Learning basada en objetos
Revisemos los archivos auxiliares que obtuvimos de la ejecución de la clasificación por random forest.
El proceso de clasificación usualmente conlleva una serie de iteraciones que implican selección de variables más importantes, búsqueda de más/mejores datos de entrenamiento y validación
Validación
- Se usan datos independientes para validar las clasificaciones
- Se construye una matriz de confusión que permite visualizar los errores por clase en los elementos que están fuera de la diagonal
- Se estiman varias medidas relacionadas a la precisión, ej.: overall accuracy y kappa
Distintas opciones: 1. Generar un nuevo set de puntos y etiquetarlos 2. Separar el set de puntos etiquetados en train y test de antemano
Opción 1
Generar un set de validación de al menos 50 segmentos. Una vez creado el vector de segmentos con etiquetas, testing, convertirlo a formato raster y ejecutar r.kappa.
r.kappa necesita mapas raster como input, por lo tanto necesitamos transformar los segmentos de validación a formato raster usando la columna class como fuente de valores para los pixeles.
Opción 2
Alternativamente, podemos separar el set de puntos etiquetados en train y test. Usemos la extensión v.divide.training_validation creada por mundialis.
Ejecutar nuevamente la clasificación usando sólo el vector training. Recordar convertir a csv.
Finalmente, convertimos a raster el vector testing utilizando la columna class como valor para los pixeles y ejecutamos r.kappa.