jueves, 21 de julio de 2011

Reportando zonas propensas a inundaciones en OpenStreetMap...

Hola Maperos!!!

Bueno, lo que pensaba sería un sencillo procedimiento para reportar las áreas inundadas (ver posts anteriores) en OpenStreetMap (OSM) no resultó tan sencillo. Decidí utilizar el script ogr2osm para convertir un archivo shapefile en un archivo con extensión .osm que pueda ser leido por JOSM, uno de los tantos clientes para manipular datos en OSM, y luego fácilmente reportarlo a sus repositorios. El problema fue que los shapefiles originados desde la clasificación de los rasters presentaban siluetas 'pixeladas' (ver figura) y el número total de nodos resultaba exageradamente alto.

Trabajando solo con la selección de los 10 más grandes cuerpos de agua, el número total era superior a los 37000 nodos. Generar el archivo .osm tomaba muchísimo tiempo y se obtenian archivos de tamaño superior a los 190 Mb. Además de que no resultaría para nada práctico subir un archivo de estas características, JOSM (o mejor, mi pobre laptop) ni siquiera era capaz de abrir un archivo de estas dimensiones. La mejor solución para este caso era usar técnicas de suavizado o generalización. La idea es eliminar un porcentaje de los nodos, sacrificando el detalle pero favoreciendo la manipulación, sin perder la información más relevante del mapa.

El algoritmo Douglas-Peucker (DP) es, tal véz, la técnica más utilizada para simplificar mapas. QGIS posee su propia implementación ([Vector -> Geometry tools -> Simplify geometries]) pero la verdad no me funcionó para nada. GRASS implementa este algoritmo, y otros, con el comando v.generalize. A travéz del plugin de GRASS para QGIS se puede tener acceso a toda su funcionalidad desde esta interfaz. Entre las opciones más interesantes esta Snake method for line smoothing (ver figura) que da una versión mucho más estilizada del mapa pero aún con gran cantidad de nodos. La version de DP en GRASS tampoco ofreció una gran mejora, solo pudo eliminar un 33% de los nodos.

La solución final fue MapSharper.org. Este es un servicio web que permite aplicar los algoritmos Douglas-Pecker y Visvalingam-Whyatt online. Aquí se puede acceder a un demo de la aplicación en Internet. Usando MapSharper la simplificación del mapa resultó trivial. El siguiente video ilustra todo el proceso:


Una vez con la version .osm del shapefile es muy simple abrirlo desde JOSM, asignarle la etiqueta flood_prone=yes y publicarlo en la base de datos de OpenStreetMap para que este a disposición de la comunidad. Se debe aclarar que antes de convertir el mapa a .osm se debe reproyectar al sistema de coordenadas WGS 84 (Grados, minutos y segundos) usado por OSM y la mayoria de aplicaciones de SIG en la web.

lunes, 18 de julio de 2011

Análisis de Imagenes Landsat TM para la Identificación de Áreas Propensas a Inundaciones.

Hola maperos!!!

Este post es para resumir el conjunto de entradas dedicadas al análisis de imagenes Landsat TM para identificar áreas propensas a inundaciones usando Software Libre. Como caso de estudio se utilizó el área del bajo cauce del rio Magdalena en el norte de Colombia durante las inundaciones del primer semestre del 2011. Aquí están los enlaces a cada uno de los tutoriales:


Happy Mapping!!!

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.