viernes, 8 de julio de 2011

Identificando cuerpos de agua y calculando estadísticas de imagenes Landsat clasificadas con SAGA GIS y QGIS

Hola maperos!!!

Llegando al final del conjunto de entradas dedicadas al análisis de imagenes Landsat para identificar áreas propensas a inundaciones en esta oportunidad veremos como tratar las imagenes clasificadas que obtuvimos en el post anterior, extraer los cuerpos de agua de cada imagen y restar uno del otro para poder cuantificar la extensión de las zonas afectadas.

La imagen clasificada que arroja SAGA GIS no es más que un raster donde cada pixel tiene un valor que representa una determinada clase (de 0 a 3 en nuestro ejemplo). Identificar a que clase pertenece cada valor será un trabajo del análista y dependerá de las características de la imagen. Por ejemplo, en algunos casos la nubosidad puede hacer parte de la clasificación o, en otros, diferentes tipos de vegetación o suelos pueden interpretarse como distintas clases por el clasificador. Sin embargo, para nuestro caso nuestro mayor interés será hacia la correcta identificación de los cuerpos de agua, el resto de las clases pueden ser omitidas.

Usando el diálogo de las propiedades de la imagen vamos a asignar un color más adecuado a cada clase en la toma clasificada para el 2011 aunque al final se extraerá la clase que se corresponde a los cuerpos de agua, que para este caso no solo incluyen rios y lagos sino también las zonas afectadas por las inundaciones. Para hacer esto usamos el comando [Raster -> Raster Calculator] en QGIS para hacer una simple operación que descarte las demás clases dejando solo la de nuestro interés. La operación:

Raster@banda = valor

hará que a todos los pixeles que correspondan al valor especificado se les asigne 1 y a todos los demás se les cambie a un valor de 0. En nuestro raster, los cuerpos de agua están representados por un valor de 2 por lo que la operación:

RioClass2011@1 = 2


será suficiente para aislar los pixeles que representan los cuerpos de agua en un nuevo raster. Solo para facilitar la visualización editaremos las propiedades de la nueva imagen para dar transparencia a los pixeles que no nos interesan y asignar un color más apropiado a los cuerpos de agua. El siguiente video ilustra lo dicho:


El proceso se repite para la imagen del 2010. Ahora contamos con dos imagenes cuyos pixeles representan con un 1 las áreas cubiertas por agua y con un 0 las que no. Usando una simple sustracción en el Raster Calculator podremos obtener aquellos pixeles que en el 2011 están cubiertos por agua (1) y no lo estaban en el 2010 (0), con lo que podriamos sugerir que se tratan de zonas afectadas por las inundaciones. El siguiente video muestra como restar de los pixeles de la imagen del 2011, sus correspondientes en el 2010:


Ahora contamos con una imagen raster donde se puede visualizar que el impacto de las inundaciones del 2011 en Colombia a sido considerable en el área de estudio, pero para cuantificar de manera precisa los alcances de las inundaciones debemos trasformar el raster en un imagen vectorial. Para hacerlo, usamos el comando [Raster -> Polygonize] en QGIS (este proceso puede tomar varios minutos). Esto nos permitirá medir el área de los polígonos resultantes usando las propiedades de la tabla de atributos que se genera junto con el archivo de polígonos (Shapefile). Esta tabla introduce un nuevo registro por cada polígono en el archivo al cual se le pueden asociar una serie de características (nombre, área, perimetro, etc.).

Para nuestro caso, primero seleccionamos solo los polígonos con valor de 1 que vienen desde el raster y al resultado adicionamos una nueva columna que calcula el área de cada polígono en Hectáreas (ha). Hay que aclarar que el Sistema de Referencia de Coordenadas original de las imagenes de Landsat (UTM Zone 18N) facilita el proceso ya que los valores calculados estarán en metros cuadrados, es por eso que solo debemos dividir por 10.000 para obtener el valor en hectáreas. Luego usamos el plugin 'Statist' para calcular las estadísticas del campo adicionado y obtener el total de hectáreas afectadas. El siguiente video muestra estos ultimos pasos:


De los resultados podemos apreciar que en total hubo alrededor de 37.241,64 ha de zonas afectadas por las inundaciones en la región y, aunque hay que hacer algunas consideraciones y limitaciones al estudio, sin duda es una cantidad preocupante. El cuerpo de agua más extenso alcanzó las 3.500 ha!!! Imagínense un solo charco de 3.500 hectáreas que antes no estaba ahí. Comparando la mayor concentración de zonas afectadas con la imagen del 2010 podemos apreciar que buena parte se trata de áreas cultivables.

Limitaciones del estudio.

Es importante también aclarar que estos valores se deben tomar con cuidado. Existen algunas limitaciones en este análisis que se deben discutir. Primero, las dos imágenes Landsat corresponden a diferentes meses del año, puede ser que no sea del todo estricto comparar los meses de Enero y Marzo. Puede que sea normal que ciertas áreas se inunden a finales de este ultimo mes. Aqui es muy importante contrastar los resultados finales con experiencia y conocimiento local para evaluar si la imagen del 2010 refleja bien el comportamiento 'natural' de los cuerpos de agua.

Otro inconveniente se da a la hora de la clasificación. La nubosidad introduce grandes problemas en esta etapa. Muchas sombras de las nubes son confundidas con cuerpos de agua por el clasificador y estos pueden impactar en el número de hectáreas afectadas. Posibles soluciones son la eliminación manual de estos 'falsos' cuerpos de agua o el uso de mascaras para eliminar las nubes y sus sombras antes de la clasificación. La conversión desde una imagen raster a sus correspondientes poligonos también puede introducir un aumento en el tamaño de los cuerpos de agua.

Bueno, con esto terminamos las entradas dedicadas al análisis de imagenes Landsat para la detección de zonas propensas a inundaciones. En un siguiente post espero poder compartir un sencillo procedimiento para reportar estas zonas en OpenStreetMap para que queden a disposición de la comunidad.

3 comentarios:

Carlos Rivas dijo...

Estimado, me ha tomado poco más de una hora seguir todo el tutorial y puedo no más que felicitarte por la claridad de los procedimientos. Espero adecuarlo a nuevas versiones de QGIS. Un saludo desde Nicaragua!

Lic. Adalberto Lujan Páez dijo...

Exelente explicación. Ahora trataré de poner en practica lo que tu con tanta sensillez has subido.
Felicitaciones por tu aporte

jorge dijo...

Hola, muchas gracias por el tutorial. No soy geografo ni nada pero he entendido el paso a paso. Tengo una pregunta como hago el calculo cuando tengo mas de 2 imagenes que comparar? En mi caso quiero ver el cambio en el tamaño del espejo de agua de una laguna en un serie de 5 años (considerando una imagen por año)

Saludos