Delimitar una cuenca hidrológica con herramientas OpenSource en Linux

Hace ya unas semanas me propuse encontrar la forma de hacer mis trabajos que regularmente hago con software privativo como ArcGIS, con herramientas Open Source y hacer una pequeña demostración de las capacidades que tienen tales aplicaciones y de la facilidad de su uso. En esta entrada presento un método para delimitar una cuenca hidrológica (Cuenca de Metztitlán en el estado de Hidalgo, México) usando la distribución ArchLinux, con la potencia de análisis que tiene el software GRASS mediante la facilidad que proporciona QuantumGIS.

Instalando las herramientas

La instalación de las herramientas necesarias para ArchLinux ya se explicó en un post anterior y realmente es muy sencilla, así que solo me limitaré a decir que las herramientas usadas son:

• ArchLinux
• GRASS 6.4
• QGIS 1.7 (con el plugin para usar GRASS)
• Además se necesita el Modelo Digital de Elevaciones de INEGI

Descargando el Modelo Digital de Elevaciones (MDE)

El modelo digital de elevaciones que se utilizó se descargó del sitio de INEGI siguiendo la ruta Geografía > Datos de relieve > Continental > Continuo de elevaciones mexicano CEM (2.0) – descarga > Descargar donde se abrirá una página en la que se tienen que introducir las coordenadas geográficas de los límites del área de estudio donde se encuentra la cuenca por delimitar. También es posible descargar la carta a escala 1:50,000 en la que se encuentra la cuenca, si se conoce el nombre de la carta (para saber la nomenclatura de las cartas visitar este sitio), y una tercera opción es descargar todo el estado en el que está la cuenca. Si se va a utilizar la segunda o tercera opción hay que asegurarse de que la cuenca se encuentre solo en una carta o en un estado, en caso de que la cuenca atraviese dos o mas cartas se pueden descargar todas ellas, lo mismo para cuando atraviesa dos o mas estados. Sin embargo es bueno descargar solo la porción que contiene la cuenca con la primera opción de descarga, la única complicación es saber las coordenadas extremas de la misma.

Descarga del MDE del sitio de INEGI

Cuando ya se hizo la selección de la descarga hay que dar clic en Procesar selección y se abrirá una página con información de la selección que se hizo anteriormente con un archivo comprimido en formato ZIP que contiene el MDE, hay que descargarlo dando clic derecho y seleccionando Guardar enlace como…. Se extrae el comprimido en un directorio y se encontrarán varios archivos:

• bil
• blw
• hdr
• prj

El archivo BIL es el que se puede abrir con QGIS y el programa lo reconocerá.

Crear una localización y conjunto de mapas

Para manejar la información y trabajar con GRASS es necesario crear una base de datos en donde se tenga el conjunto de capas deseadas. Una base de datos de GRASS es simplemente un directorio en el que se guardan locations (localizaciones) y mapsets (conjuntos de mapas). Se pueden tener todas las bases de datos que se deseen, para usar GRASS generalmente se crea una llamada grassdata en la carpeta personal: home/usuario.

En este caso se creará la base de datos con el plugin de QGIS-GRASS el cual permite usar los módulos de GRASS desde el ambiente amigable de QGIS. Para activar el plugin se abre QGIS y se va al menú Complementos > Administrar complementos… en donde se tiene que marcar el que se llama GRASS, yo personalmente activo todos.

Activando el plugin de GRASS en QGIS

Cuando se tenga activado el complemento se puede seguir con la creación de la base de datos. Dar clic en el icono Añadir capa ráster en la barra de menú de QGIS y buscar el ráster que se descargó de INEGI y seleccionar el archivo en formato BIL y abrirlo. Una vez cargado, para crear la base de datos selecciona el menú Complementos > GRASS > Nuevo directorio de mapas y se abrirá una ventana que pedirá crear o seleccionar una carpeta para la base de datos.

Creando una nueva base de datos

Se pide elegir el directorio para crear la base de datos, como se mencionó anteriormente se elige crear un directorio llamado grassdata en la carpeta de usuario.

Elegir directorio para la base de datos

Se crea una nueva location para ello simplemente se le da el nombre deseado, en este caso se le dió el nombre del estado en el que se encuentra la mayor parte de la cuenca.

Crear nueva localización

Se elije el sistema de coordenadas de referencia,las cuales pueden ser coordenadas geográficas o elegir un sistema de proyección, en la imagen se muestra que el seleccionado es el sistema de proyección UTM con elipsoide y datum WGS84 y Zona 14N, la elección depende de los objetivos del usuario. Para ayudar a elegir la zona UTM se puede hacer uso de GoogleEarth, marcando el menú Ver > Cuadrícula y luego en Herramientas > Opciones y marcar Universal Transverse of Mercator y se mostrarán las zonas UTM en el globo. Elegir la que corresponda.

Selección de la proyección

A continuación se elige la región predeterminada para trabajar, en este caso elegí la que corresponde para México de la lista.

Región predeterminada de GRASS

Finalmente aparece una ventana en la que se pide el nombre del mapset o conjunto de mapas (layers) del usuario, introducir el nombre deseado (p. ej. el nombre del usuario) y luego se pide confirmación de los datos. Al aceptar se establecerá automáticamente la nueva localización y conjunto de datos por defecto para QGIS.

Con estos sencillos pasos se tiene una nueva base de datos creada, al principio parece difícil pero es cuestión de acostumbrarse a la manera en que GRASS organiza las cosas. Usando QGIS el proceso de creación de la base de datos ha resultado un poco más fácil que si se hubiera utilizado solo GRASS.

Importar ráster a la base de datos

Ya que se ha creado la base de datos, se puede notar que está vacía puesto que la capa ráster en formato BIL agregada al principio, está en otro directorio fuera de la nueva base de datos. Para importar esta capa o layer a la base de datos se utiliza la opción: Complementos > GRASS > Abrir herramientas de GRASS y en la ventana que aparece escribir y seleccionar el módulo r.in.gdal.

Importar ráster GDAL

En la ventana del módulo que aparece seleccionar el ráster a importar y darle un nombre de salida. También puede consultarse el manual en inglés en la pestaña Manual. Para procesar la imagen dar clic en Ejecutar, y al terminar el proceso en Ver salida para que se agregue la capa al QGIS automáticamente. La capa importada se agrega al mapset de la base de datos por defecto, es decir la última creada.

Ráster importado y agregado de manera automática

Delimitación de la cuenca

Para la delimitación de la cuenca se utiliza el módulo r.watershed, abriéndolo desde Complementos > GRASS > Abrir herramientas de GRASS y seleccionándolo en la lista de módulos, aparecerá una ventana en la que se introducen los datos para necesarios.

Módulo r.watershed ejecutándose

Se elije como entrada el ráster de elevaciones importado llamado «mdemetztitlan» el cual aparece mientras esté marcada la casilla para ser visto en la sección Capas, todo el análisis se basará en este ráster. Un dato importante es el tamaño mínimo para cada cuenca el cual se tiene que introducir en número de celdas, para éste ejemplo se colocó un tamaño mínimo de 111 celdas, haciendo el cálculo como sigue:

Resolución del ráster de INEGI (una celda): 30m x 30m = 900 m² entonces para una hectárea que equivale a 10,000 m² entre 900 m² = 11.11 celdas / ha. Entonces si se desean delimitar cuencas con un tamaño mínimo de 10 ha se tendrá que introducir un valor de 111.11 celdas.

Se pide introducir los nombres de las capas de salidas que para este ejemplo son:

acumulacion
Mapa de salida: El valor absoluto de cada elemento en esta capa de salida es la cantidad de flujo superficial que atraviesa la celda. Este valor será el número de celdas terreno arriba, más una, si no se especifica el mapa de flujo superficial. Si se da el mapa de flujo superficial, el valor estará en las mismas unidades que flujo superficial. Los números negativos indican que las celdas posiblemente tienen escurrimiento superficial desde fuera de la región geográfica actual. Por lo tanto, en celdas con valores negativos no se puede calcular con precisión el escurrimiento superficial y la producción de sedimentos.

drenaje
Mapa de salida: la dirección de drenaje. Proporciona el «aspecto» de cada celda medido en sentido contrario a las manecillas del reloj desde el Este. Multiplicando los valores positivos por 45 dará la dirección en grados que el escurrimiento superficial viajará desde la celda. El valor 0 (cero) indica que la celda es una zona de depresión (definida por el mapa de entrada de depresiones). Los valores negativos indican que el escurrimiento superficial se sale de los límites de la región geográfica actual. El valor absoluto de estas celdas negativas indica la dirección del flujo.

corrientes
Mapa de salida: segmentos de la corriente. Los valores corresponden a los valores de las cuencas hidrográficas. Pueden ser vectorizados después de un adelgazamiento (r.thin) con r.to.vect.

cuencas
Mapa de salida: da una etiqueta única para cada cuenca hidrográfica. A cada cuenca se le asignará un único número par entero y positivo. Áreas a lo largo de los bordes no pueden ser lo suficientemente grandes como para crear una cuenca hidrográfica exterior. Valores de 0 indican que la celda no es parte de una cuenca hidrográfica completa en la región geográfica actual.

Dar clic en Ejecutar y al terminar el proceso dar clic en el botón Ver salida para agregar automáticamente el resultado a QGIS. De manera adicional se muestra el comando que se ejecuta en el shell de GRASS, el cual se podría ejecutar si únicamente se usara ésta aplicación, en QGIS se pone con una finalidad informativa:

r.watershed elevation=mdemetztitlan@Eduardo threshold=111 accumulation=acumulacion drainage=drenaje stream=corrientes basin=cuencas

Resultando:


SECTION 1a (of 5): Initiating Memory.
SECTION 1b (of 5): Determining Offmap Flow.
SECTION 2: A * Search.
SECTION 3: Accumulating Surface Flow with SFD.
SECTION 4: Watershed determination.
SECTION 5: Closing Maps.
Finalizado correctamente

En la corrida anterior se usó un valor mínimo de cuencas muy pequeño y por ello como resultado se obtuvo una capa con muchas cuencas pequeñas, para delimitar la cuenca «principal» es conveniente entonces incrementar este valor y volver a ejecutar el módulo. Ahora para delimitar la cuenca general se usa un tamaño mínimo de cuencas de 10,000 ha ejecutando el módulo con un área de 111111 celdas se obtiene un resultado como el mostrado en la imagen:

Delimitación de la cuenca de Metztitlán

Clic en el botón Ejecutar y cuando termine en el botón Ver salida para agregar automáticamente el resultado a QGIS. El código de la instrucción ejecutada en GRASS ahora es:

r.watershed elevation=mdemetztitlan@Eduardo threshold=111111 accumulation=acumulacion2 drainage=drenaje2 stream=corrientes2 basin=cuenca2

El resultado que se despliega es:

SECTION 1a (of 5): Initiating Memory.
SECTION 1b (of 5): Determining Offmap Flow.
SECTION 2: A * Search.
SECTION 3: Accumulating Surface Flow with SFD.
SECTION 4: Watershed determination.
SECTION 5: Closing Maps.
Finalizado correctamente

De esta forma se obtuvieron dos resultados, por una parte, al ejecutar el módulo con un tamaño mínimo de cuencas de 10 ha, se obtuvieron cuencas con una superficie pequeña y una red de corrientes detallada que más adelante servirá para generar una capa vectorial de escurrimientos superficiales, además de la capa de drenaje y acumulacion. Y por otra parte se obtuvo una delimitación general de la cuenca usando un tamaño mínimo de 10,000 ha, en este caso se nombraron las salidas con los mismos nombres que la anterior pero agregando un «2» al final quedando acumulacion2, corrientes2, cuenca2 y drenaje2.

Ahora para manejar de una forma más práctica la delimitación de la cuenca, se convierte la capa llamada cuenca2, de ráster a vectorial. Para hacer la conversión hay que abrir las herramientas de GRASS y hacer clic en la pestaña Lista de módulos en la que hay que seleccionar r.to.vect.area el cual sirve para convertir un archivo ráster a vectorial (polígonos).

Conversión de ráster a vectorial

El código de la instrucción es:

r.to.vect input=cuenca2@Eduardo output=cuencametztitlan feature=area -s

Se obtiene la siguiente información al ejecutar el módulo, notar que salen un par de advertencias pero nada de cuidado por lo pronto:

!Default driver / database set to:
!driver: dbf
!database: $GISDBASE/$LOCATION_NAME/$MAPSET/dbf/

Extracting areas...
Building topology for vector map ...
Registering primitives...
2 primitives registered
4192 vertices registered
Building areas...
1 areas built
1 isles built
Attaching islands...
Attaching centroids...
Number of nodes: 2
Number of primitives: 2
Number of points: 0
Number of lines: 0
Number of boundaries: 1
Number of centroids: 1
Number of areas: 1
Number of isles: 1
r.to.vect complete.
Finalizado correctamente

Al terminar la ejecución de la instrucción, se tiene la delimitación de la cuenca en vectorial. Resulta más práctico manejar áreas en formato vectorial como polígono que en vez de un ráster, además en mi opinión se le puede dar más presentación a un polígono.

Para terminar hay que hacer un sombreado de las elevaciones, conocido como hillshade que tiene fines estéticos. Para ello usar el módulo llamado r.shaded.relief que requiere lo siguiente: el modelo de elevaciones, en este caso mdemetztitlan; el nombre de la capa de salida, yo lo llamé sombreado; la altitud del sol en grados por arriba del horizonte, el azimut del sol en grados al Este del Norte y la exageración, los cuales dejé lo que pone por defecto. Dar en Ejecutar y al terminar en el botón Ver salida para agregar automáticamente la capa resultante a QGIS.

Creando el sombreado

Listo ahora agregar un poco de transparencia para obtener el resultado final como se muestra en la imagen.

Delimitación de la Cuenca de Metztitlán

Generar corrientes

Para generar corrientes superficiales se utiliza la capa corrientes generada al principio del análisis puesto que la red es mucho más detallada que corrientes2, esta capa es ráster y hay que convertirla a vectorial. Antes de ejecutar la conversión es necesario procesarla por medio de la herramienta r.thin como se muestra en la imagen:

Adelgazando la capa ráster con r.thin

La capa de salida se llama rios misma que se convierte a vectorial con el módulo r.to.vect.line de GRASS. El procedimiento únicamente pide el nombre de la capa de salida que en este ejemplo de llama escurrimientos.

Generando escurrimientos superficiales

Cuando termine, se tendrán escurrimientos generados con un área de captación de 10 ha, para obtener menor o mayor detalle se puede repetir el proceso, cambiando el área mínima de la cuenca por generar.

Corrientes generadas

Finalmente resta guardar el proyecto llendo al menú Archivo > Guardar proyecto como … y guardarlo en el directorio y con el nombre que se desee.

Conclusiones

Con este ejemplo pretendí hacer una pequeña demostración de la eficiencia y calidad que tienen las herramientas Open Source que existen para el manejo de sistemas de información geográfica, con las que se pueden obtener muy buenos resultados.

En mi opinión, en este pequeño «artículo» se pudo notar la potencia que tiene la aplicación GRASS por medio de sus módulos que hacen una gran cantidad de tareas y que además las hacen bien, sin embargo el uso de este paquete podría resultar complicado para usuarios principiantes que vienen de usar software privativo, por ello una buena opción es manejar los módulos de GRASS por medio de QGIS para facilitar la tarea.

Por último resta decir como profesionistas no podemos quejarnos de las herramientas libres que se ofrecen puesto que resultan bastante buenas, incluso algunas son mucho mejores que las soluciones de software privativo además, en mi opinión personal, académicamente es más conveniente usar aplicaciones en las que se pueda tener acceso al código fuente porque se puede saber exactamente que algoritmos de cálculo o modelación se usan, aparte de la ventaja de no tener que pagar costosas licencias por usar determinados paquetes de software.

Delimitación final de la Cuenca de Metztitlán

¡Saludos y espero les sirva!

Fuentes:

http://earth.unibuc.ro/tutoriale/bazin-hidrografic-prin-QGIS-GRASS/
Manual de GRASS

NOTA: La cuenca no está completa pero sirve como ejemplo para la realización de la delimitación.