En los últimos meses que estuve inmerso aprendiendo Python (razón principal que explica mi ausencia en el blog), me tope con una interesante herramienta que incluye una versión como un plugin del QGIS. Como su desarrollo está basado enteramente en Python, me permitió explorarlo un poco más de lo habitual y quisiera poder compartirles dicha experiencia.
Me estoy refiriendo al denominado "Relief Visualization Toolbox (RVT)", el cual según sus creadores, tiene como objetivo ayudar a los científicos a visualizar conjuntos de datos de modelos de elevación ráster. Señalan que la configuración predeterminada asume el trabajo con modelos digitales de elevación de alta resolución, derivados de misiones de escaneo láser aerotransportadas (LIDAR). Es importante mencionar que esta herramienta fue posible gracias al trabajo de un grupo de investigadores del Institute of Anthropological and Spatial Studies.
Librería de Python rvt-py
Lo primero que debemos hacer es tener instalado la librería de Python rvt-py, en mi caso lo hice siguiendo las indicaciones del pypi (Python Package Index). Desde nuestra consola solo debemos tipiar la siguiente instrucción.
pip install rvt-py
Es importante que antes de hacer eso, verificar que tengamos instalado las librería de numpy, scipy y gdal.
Vamos a describir los pasos a seguir para obtener un grupo de visualizaciones que la herramientas es capaz de generar partiendo de un modelo de elevación digital (DEM), todo ello apoyado de la documentación disponible de la herramienta (1). Es importante mencionar que esta librería de Python contiene dos módulos principales: Modulo vis y Modulo default.
Paso 1: Inicio. En nuestro editor de código Python vamos a configurar los requerimientos importanto las herramientas necesarias.
Paso 2: Cargando nuestro DEM. Aquí vamos a cargar nuestro archivo DEM en formato GeoTiff a una matriz numpy, para ello usaremos el módulo rvt.default (que se basa en la biblioteca gdal de python). También podríamos usar rasterio, gdal o cualquier otra biblioteca de Python que tenga esa capacidad. Para nuestra prueba vamos a emplear un DEM de una zona desértica del Municipio de Mexicali, dentro del Estado Mexicano de Baja California. Lo importante es que tiene un tamaño de píxel de 0.5 m.
El módulo default tiene la función get_raster_arr () que lee el raster según la ruta indicada y devuelve un diccionario con las claves "array", "resolución" y "no_data". La clave "array" es un numpy array del ráster, la "resolución" es una tupla del tamaño de píxel donde el primer elemento es el tamaño de píxel en la dirección 'x' y el segundo en la dirección 'y', "no_data" es el valor que representa noData en el raster (array).
Paso 3: Generando las visualizaciones: Ahora que ya tenemos a nuestro DEM el cual representa el parámetro de entrada, podemos ejecutar las funciones de visualización disponibles con la librería, entre ellas tenemos a: Pendiente, Hillshade (sombreado), Hillshade de múltiples direcciones, Modelo de relieve local simple, Factor de vista del cielo, Factor de vista del cielo anisotrópico, Apertura - Positivo, Apertura-Negativo, Iluminación del cielo y Dominio local. Cada uno de ellos toman algunos parámetros que son específicos para cada visualización y parámetros comunes, los cuales iremos detallando.
1. Pendiente - Aspecto.
Para calcular la pendiente y el aspecto, debemos usar la función rvt.vis.slope_aspect (). Los parámetros son: dem, resolution_x, resolution_y, output_units (puede ser: percent, degree, radianes), ve_factor (factor de exageración vertical) el valor por defecto es 1, pero si presenta un relieve plano se recomienda usar > 1 , no_data (valor que representa no_data, todos los píxeles con este valor son cambiados a np.nan), fill_no_data (bool, si es 'True', todos los píxeles no_data se reemplazan con la media de píxeles vecinos), keep_original_no_data (Si es True, cambia todos los píxeles de salida a np.nan donde el DEM presenta no_data).
Para guardar la visualización usaremos rvt.default.save_raster (). Esta función presenta los parámetros: src_raster_path (ruta dem para copiar geodatos), out_raster_path (ruta del nuevo archivo, en donde se guardará nuestro ráster generado), out_raster_arr (matriz con datos ráster), no_data (cómo los 'no datos' están almacenados, todas las visualizaciones no_data se almacenan como np.nan), e_type (GDALDataType, por ejemplo, 6 es para Float32 y 1 es para UInt8).
Lo podemos visualizar en QGIS.
2. Hillshade (Sombreado)
Para enteder un poco mejor sobre el sombreado debemos tener claro que representan azimut y elevación. Según (2), son las dos coordenadas que definen la posición de un cuerpo celeste (sol, luna) en el cielo visto desde una ubicación particular en un momento particular.
 |
| Representación del azimut y elevación del sol. (Fuente: photopills.com) |
Para calcular el sombreado usaremos la función rvt.vis.hillshade (). Los parámetros nuevos son: sun_azimuth (315 es el valor por defecto, siendo el rango 0-360), sun_elevation (35 es el valor por defecto, siendo el rango 0-75) se recomienda usar valores pequeños (5,10) para terrenos planos y valores altos (45) para terrenos empinados. La función genera un "numpy array 2D" de sombreado.
Ahora solo nos queda guardar nuestros resultados para visualizarlo en nuestro QGIS.
3. Sombreado Multidireccional
A diferencia de un sombreado predeterminado, el sombreado multidireccional presenta una vista incomparable de las montañas, mesetas, valles y cañones del mundo mediante el uso de un algoritmo que calcula el sombreado desde seis direcciones diferentes en lugar de la dirección única utilizada en un sombreado predeterminado. El resultado es una visualización impresionante tanto en pendientes altas como en áreas inexpresivas (3). Mejora el equilibrio entre las áreas sobreexpuestas y las áreas de sombra del mapa. La salida es adecuada para usarla como telón de fondo en relieve para los mapas topográficos, de suelo, hidrológicos, de cobertura de suelo y otros mapas temáticos en los que los datos se realzarán con la topografía (4).
Para calcular el sombreado de múltiples direcciones, use la función rvt.vis.multi_hillshade (). El parámetro nuevo a ingresar es: nr_directions (16 es el valor por defecto en un rango de 4-360) recomiendan no ingresar valores mayores a 16 debido a la alta autocorrelación. La función produce un numpy array 3D (donde la primera dimensión representa cada dirección (nr_directions), por ejemplo, arr[0] es la primera dirección).
Tener en cuenta que al guardar la matriz de sombreado de múltiples direcciones, cada dirección (azimut solar) se guardará en una banda.Si revisamos en QGIS las propiedades de nuestra capa ráster veremos las 16 bandas generadas.
4. Modelo de Relieve Local Simple (SLRM)
El LRM representa las diferencias de elevación locales a pequeña escala después de eliminar las formas del paisaje a gran escala de los datos. El LRM mejora en gran medida la visibilidad de las características topográficas superficiales a pequeña escala, independientemente del ángulo de iluminación, y permite medir directamente sus elevaciones relativas y sus volúmenes (5).
El modelado de relieve local elimina los elementos morfológicos a gran escala (colinas, valles ...) de los datos, por lo que solo quedan características a pequeña escala (por ejemplo, arqueología). Esta versión de la herramienta utiliza un proceso simplificado: la tendencia se calcula mediante un filtro de media simple y un modelo eliminado de tendencia se produce directamente restando el modelo filtrado del original (6).
Para calcular el modelo de relieve local simple, use la función rvt.vis.slrm (). El nuevo parámetro a considerar es radius_cell (20 px es el valor por defecto en un rango de 5 -50) tener en cuenta que el radio debe ser un poco más de la mitad de los elementos que le interesan. La función devuelve una gran variedad de modelos de relieve locales simples.
5. Modelo de Relieve Multi-Escala
El Modelo de relieve multi-escala (MSRM - Multi-Scale Relief Model), es un nuevo algoritmo para la interpretación visual de accidentes geográficos utilizando DSM (Modelos de Superficie Digital).La importancia de este nuevo método radica en su capacidad para extraer la morfología del relieve de los DSM de alta y baja resolución, independientemente de la forma o escala del relieve en estudio (7).
Para calcular el modelo de relieve de múltiples escalas, utilice la función rvt.vis.msrm (). Los parámetros nuevos a considerar son: feature_min, feature_max y scaling_factor. En relación al factor de escala, si es mayor que 1, proporciona un rango más amplio de valores de MSRM (aumenta el contraste y la visibilidad), pero podría resultar en una pérdida de sensibilidad para los elementos de tamaño intermedio. La función devuelve una gran variedad de modelos de relieve de múltiples escalas.
En esta primera parte se presentaron 5 tipo de visualizaciones que podemos generar empleando la libreria de Python rvt-py. Existen un grupo de otras opciones que serán presentado en una segunda parte, del mismo modo se mostrará la posibilidad que tenemos de hacer combinaciones.
Plugin de QGIS Relief visualization toolbox (RVT)
Ahora que hemos revisamos con cierto detalle las opciones para generar visualizaciones, luego de instalarlo en nuestro QGIS va ser más fácil de entenderlo. En realidad es muy sencillo su uso, solo debemos contar con un DEM en nuestro panel de capas, luego activamos la herramienta y seleccionamos aquellas visualizaciones que deseamos contar, podemos dejar con las valores predeterminadas para los parámetros que lo requieran, o cambiarlos ahora que conocemos los rangos disponibles para la mayoría de ellos.
Recuerda que puedes seleccionar una carpeta de salida específica, para ello debes quitar el check de la casilla 'save to raster location'. Si seleccionaste todas la visualizaciones y luego hiciste click en el botón 'Start', luego de unos minutos (en este caso demoró un poco más de 5 minutos) tendras un grupo de archivos ráster que podrás verlo en el QGIS.
Como seguro lo vieron la herramienta tiene otra pestaña denominada 'Blender' con el que podemos hacer combinaciones, pero de eso hablaremos en la segunda parte.
Por ahora se comparte un repositorio (https://github.com/ccarbajal16/python_rvt_test) que se creó para ir avanzando en las pruebas con esta herramienta, ahí encontrarán el código en Python asimismo el avance de un Notebook de Jupyter. Se aclara que el DEM original que se muestra aquí fue recortado con fines prácticos.
Referencias:
- Relief Visualization Toolbox in Python (https://rvt-py.readthedocs.io/en/latest/)
- Understanding-azimuth-and-elevation (https://www.photopills.com/articles/understanding-azimuth-and-elevation)
- Multi-Directional Hillshade Makes Your Maps Pop (https://www.esri.com/about/newsroom/arcuser/multi-directional-hillshade-makes-your-maps-pop/)
- Función Sombreado (https://pro.arcgis.com/es/pro-app/latest/help/analysis/raster-functions/hillshade-function.htm)
- Using LIDAR-derived Local Relief Models (LRM) as a new tool for archaeological prospection (https://www.jstor.org/stable/j.ctt6wp79m.30)
- RVT manual, version 2.2.1 (https://www.zrc-sazu.si/sites/default/files/rvt_2.2.1_0.pdf
- Multi-scale relief model (MSRM): a new algorithm for the visualization of subtle topographic change of variable size in digital elevation models.(https://onlinelibrary.wiley.com/doi/full/10.1002/esp.4317)
No hay comentarios.:
Publicar un comentario