Recomendaciones para el pre-procesamiento de modelos digitales de elevación en análisis de cuencas

Antes de comenzar, quisiera indicarte que este artículo sólo muestra recomendaciones basadas en las experiencias de manejo de HidroCalc 2 en los procesos de delimitación de cuencas y obtención de sus parámetros morfológicos. Debes tomarte la información como una ayuda y no como un manual a rajatabla. Si no estás de acuerdo con la información aquí vertida o consideras que se puede mejorar en algunos aspectos, siéntete libre de manifestar tu opinión en la caja de comentarios o a través del correo unu.ingenieros@gmail.com, si nos ponemos de acuerdo, mejoraremos el artículo y colocaremos tu nombre en los créditos.

Como haz de haber leído por ahí, un modelo digital de elevación (MDE o DEM por sus siglas en inglés), puede provenir de diferentes fuentes: levantamientos topográficos, vuelos fotogramétricos o satélites. Estos últimos son los que comúnmente utilizamos para la delimitación de cuencas hidrográficas y los podemos obtener de diferentes agencias espaciales a través de las páginas web de los organismos gubernamentales de sus respectivos países, no entraré en más detalles ajenos a nuestro propósito.

A nosotros, estos MDE pueden llegar en diferentes formatos y niveles de procesamiento y por lo tanto cada uno de ellos requiere de un tratamiento especial según el propósito para el que lo utilizaremos, en este caso, para delimitar una cuenca y obtener sus parámetros morfológicos. En cuanto al formato, me limitaré a los MDE de grillas cuadradas.

Antes de empezar, debo hacer algunas precisiones respecto a lo que representa un MDE y un MDT.

Respecto al MDE: es un modelo de coordenadas (x,y,z) de una superficie almacenada digitalmente y es necesario especificar qué tipo de superficie es representada: superficie de vegetación, superficie de agua subterránea, superficie del terreno, etc.

Respecto al MDT: es un modelo de coordenadas triples (x,y,z) de la superficie terrestre almacenada digitalmente, convirtiéndolo en un caso particular del MDE.

El término "superficie terrestre", puede referirse a la superficie natural y la hecha por el hombre y se limita a la capa entre la atmósfera o los cuerpos de agua y la corteza terrestre.

Hasta aquí, de lo que se trata este artículo es de cómo convertir un MDE a un modelo digital de terreno (MDT o DTM por sus siglas en inglés),.

Los procesos comúnmente aplicados son: transformación del MDE a un sistema de coordenadas proyectadas, relleno de vacíos y quemado de la red de drenaje.

Sin embargo, existen otros procesos no tan comunes, aplicados en ámbitos más cercanos a la investigación en teledetección y software para SIG o situaciones en las que se requiere mayor grado de precisión, tales como: eliminación de ruido, limpieza de características antrópicas y limpieza de bosques. Por supuesto, la información requerida y los algoritmos de procesamiento son más complejos para estas tareas que para las anteriores.

A continuación, procederé a explicar brevemente de qué se trata cada uno.

Proyección: este proceso consiste en transformar el MDE desde un sistema de coordenadas esféricas a un sistema de coordenadas proyectadas, esto es debido a que existe la necesidad de calcular longitudes y áreas en los procesos de delimitación y obtención de características morfológicas de la cuenca.
Reducción de ruido: el MDE obtenido de radar incluye una variación local de la elevación, esta variación es conocida como ruido térmico y no contiene ninguna información relevante en el análisis del terreno, por tal motivo, es necesaria su reducción.
Relleno de vacíos/datos faltantes: dada la forma de captura de información, las plataformas que portan los sistemas de radar suelen generar cierto ángulo con la superficie terrestre que según la posición en la que se encuentren, puede incidir en superficies elevadas y perder, o no capturar, la información que se encuentra en la parte posterior. Esto genera espacios vacíos en el MDE, y por esta misma razón, suelen presentarse en superficies con pendientes y elevaciones grandes. Esta pérdida de información impide el análisis completo del terreno, por ello, es importante realizar el relleno.
Limpieza de características antrópicas: para analizar el relieve en condiciones naturales es necesario eliminar del MDE la superficie modificada por el hombre, por ejemplo, represamientos u otras grandes estructuras. Un algoritmo para realizar este proceso no requiere por lo general de información adicional dado que se utiliza un análisis de tendencia local en el que se elimina la información correspondiente a elevaciones atípicas.
Limpieza de bosques: consiste en eliminar o filtrar del MDE las grandes masas forestales. A diferencia del proceso anterior, lo que aquí se necesita es la ubicación del bosque.
Quemado de la red de drenaje: este proceso consiste en forzar que el MDE presente depresiones en donde en la realidad existen cauces o ríos. Este proceso es necesario debido a la pérdida de información que pueden acarrear los procesos anteriores o a la misma precisión del MDE. En zonas planas o de bajas pendientes o en cauces poco profundos, es muy importante realizar este proceso para obtener una correcta delimitación del MDE.

De lo expuesto anteriormente, se puede deducir que cuando se trata de obtener las características morfológicas de una cuenca hidrográfica, lo más recomendable es realizar todos los procesos, en el orden que fueron presentados.

Dado que son pocas las situaciones en las que se dispone de información o herramientas para realizar el proceso 4, podemos llegar a omitirlo, siempre que el nivel de precisión del MDE y las circunstancias lo permitan, teniendo en cuenta que muchas cuencas en nuestro territorio no presentan severas características antrópicas o de hecho no las presentan y que por lo general no disponemos de MDE con la precisión necesaria para evaluar la magnitud de las características antrópicas.

El paso 5 es sumamente importante, pero nuevamente puede omitirse considerando el nivel de precisión del MDE y los requerimientos del análisis. En mi opinión, el hecho de no disponer de información suficiente no debe ser un motivo para saltarse este proceso en MDE's de resolución espacial mayor a 30 m, considerando que las masas forestales pueden ser obtenidas a partir de imágenes multiespectrales de libre distribución, tales como las que provienen de Sentinel o Landsat, cuyas resoluciones pueden llegar incluso a los 10 y 15 m, respectivamente, mediante procesos de fusión.

Los procesos restantes ya poseen herramientas para su análisis en software SIG, como QGIS o ArcMap. No es deberías omitir ninguno si buscas realizar un análisis adecuado del terreno, sin embargo, el proceso 3, en cierto grado, es menos relevante que los tres primeros.

En este artículo, mostraré las herramientas necesarias para realizar los procesos 1,2,3 y 6. Los procesos 4 y 5 los dejaré para un artículo futuro.

Para esto, utilizaré el MDE más popular, el SRTM, específicamente en su producto de 1 arco segundo (30 m), para explicar cómo realizar los procesos que deberás extrapolar a los demás.

Proyectar: Este paso puedes realizarlo en QGIS con la herramienta 'Combar (reproyectar)' de GDAL, o en ArcMap mediante la herramienta 'Project Raster' de ArcToolbox, sobra decirte que puedes escribir código en Python para este proceso.

Reducción de ruido: En QGIS se puede ejecutar la eliminación de ruido con la herramienta de SAGA 'Mesh denoise' e ingresando solamente el MDE. En ArcMap el proceso se puede hacer con la herramienta 'Filter'. Para que puedas notar los cambios, he aplicado a este MDE una visualización del tipo 'Hillshade'.

La siguiente animación muestra el pre, post y la diferencia entre ambos, luego del proceso.

La diferencia entre los MDE no es tan severa como para descartar el proceso por asumir que se ha alterado significativamente los datos, aquí puedes ver el histograma de la diferencia del pre y post:

Este proceso es más notorio cuando se realiza en un MDE de mejor resolución espacial y cuando se ejecutan un mayor número de iteraciones. Aquí te muestro el resultado de proceso de un Alos Palsar, y un GIF que obtuve de stackexchange:

Rellenar faltantes: Ahora necesitamos completar la información faltante en el MDE, para esto podemos utilizar directamente herramientas programadas en QGIS, tales como 'r.neighbors', o en ArcMap con la herramienta 'Image Analysis" y el uso de su máscara de función 'Elevation Void Fill Function' (solamente disponible en versiones superiores a la 10.3, aquí tienes un video). Si dispones de versiones anteriores, puedes utilizar la herramienta 'Focal Statistics'.

Debes tener en cuenta que este proceso puede ocasionar pérdida de precisión del MDE. De hecho al finalizar, tendrás que aplicar un filtro de selección mediante la calculadora ráster y copiar solo la información completada hacia tu MDE original si es que haz realizado muchos pasos. También hay que tener en cuenta que el proceso de eliminación de ruido suele rellenar un poco los espacios vacíos.

Es mejor recortar del MDE a tu zona de interés antes de hacer el proceso. Se debe indicar el "tamaño del vecindario de celdas", debe ser un número impar y colocar números muy grandes no da buenos resultados. La herramienta se repite tantas veces como sea necesario hasta rellenar los espacios vacíos, difícilmente se habrán llenado todos los espacios la primera vez.

En esta animación puedes ver que utilice 16 veces la herramienta (a modo de ejemplo), colocando cada vez el resultado del proceso anterior.

Quemado de la red de drenaje: Para este proceso, es necesario disponer de la red de ríos y lagos de la cuenca. Esta información se puede obtener de mapeo en campo o en su defecto, de fotografía aéreas. En Perú, se dispone de la información cartográfica a escala 1:100'000 dada por el IGN (aquí el enlace de descarga) que contiene, entre otra información, ríos y lagos. Lamentablemente, los datos no tienen elevación, por lo que no podemos sacar mucho provecho de ellos. Aun así, podemos utilizar la información de ríos y será de mucha ayuda sobre todo en zonas planas.

El quemado de la red de drenaje, es más fácil de realizar en ArcMap con las herramientas de ArcHydroTools, específicamente 'DEM reconditioning'. Sin embargo, también en QGIS se puede realizar el proceso con la herramienta de SAGA 'Burn stream network into dem' o de GRASS con 'r.carve'. La herramienta de SAGA requiere de la red de drenaje en formato de ráster, lo que complica su uso. La de GRASS es bastante más similiar a la de ArcHydroTools.

Algunas consideraciones adicionales que te pueden servir en base a mi experiencia son las siguientes:

Verifica que el MDE, tu archivo de ríos y la configuración del software para SIG se encuentren en el mismo sistema de coordenadas proyectadas.
Es probable que en algunas cuencas necesites unir más de una red de drenaje de varias cartas, para ello debes unir los archivos *.shp mediante las herramientas 'Merge'.
Si deseas añadir o eliminar ríos en tu archivo (algo muy común dada la antigüedad y la escala de las cartas IGN), no debes olvidar verificar la topología del archivo. Puedes hacerlo manualmente o utilizar las herramientas topológicas de los software para SIG, es probable que para esta verificación requieras mayores conocimientos sobre geodatabases.
El ancho de los cauces puede ser menor en la realidad, sin embargo, colocar el verdadero ancho podría no dar buenos resultados en este algoritmo. Debes considerar el ancho de las celdas y permitir que se pueda modificar en al menos 2. Esto es, a pesar de que el cauce del ejemplo tiene un ancho promedio de 8 metros, el MDE tiene un tamaño de celda de 30 metros, por lo tanto, o se cambia dicho tamaño, o se coloca un ancho de 60 metros en el algoritmo.

Finalmente, quiero compartir algunos procesos que suelo hacer a los MDE que utilizo con mayor frecuencia.

Descarga SRTM - Mundial