Uso de datos Lidar GEDI con QGIS apoyados con Python

 

En esta oportunidad motivado por poner en práctica una capacitación brindada por ARSET, les quiero mostrar el procedimiento que nos permite contar con los datos Lidar GEDI (Global Ecosystem Dynamics Investigation), instrumento que está a bordo de la Estación Espacial Internacional (ISS) y su misión tiene como objetivo caracterizar la estructura y la dinámica del ecosistema para permitir una cuantificación y una comprensión radicalmente mejorada del ciclo del carbono y la biodiversidad de la Tierra (1). 

Como lo mencionan en https://gedi.umd.edu/, GEDI proporcionará respuestas sobre cómo la deforestación ha contribuido a las concentraciones de CO2 atmosférico, cuánto carbono absorverá los bosques en el futuro y cómo la degradación del hábitat afectará la biodiversidad global. Todo ello gracias a que GEDI logra observar casi todos los bosques tropicales y templados utilizando un altímetro láser autónomo de la ISS. 

El instrumento GEDI es un sistema láser de detección de luz y rango (lidar) de clase geodésica compuesto por 3 láseres que producen 8 pistas paralelas de observaciones. Cada láser dispara 242 veces por segundo e ilumina un punto de 25 m (una huella) en la superficie sobre la que se mide la estructura 3D. Cada huella está separada por 60 m a lo largo de la vía, con una distancia entre vías de aproximadamente 600 m entre cada una de las 8 vías. Es importante mencionar que recopila datos a nivel mundial entre las latitudes 51.6 ° N y 51.6 ° S. Una animación de cómo funciona lo encuentran aquí.

Patrón de muestreo terrestre de GEDI
(Fuente: https://gedi.umd.edu/instrument/specifications/)

Productos de GEDI:

En resumen con GEDI podremos tener información sobre la altura del dosel del bosque, la estructura vertical del dosel y la elevación de la superficie. Para ello están disponibles digamos tres tipos de productos, todos ellos tienen una resolución espacial de 25 m y son brindados en formato HDF5 (*.h5).

• GEDI Level 1B: Producto de formas de ondas geolocalizadas (GEDI01_B). Contiene 85 capas para cada uno de los ocho haces. Deben considerar que tiene un tamaño de archivo aproximado de 7GB. Para mayor información pueden consultar aquí, también existe una guía de usuario del producto en pdf.

Formas de ondas geolocalizadas
(Fuente: Elaboración Propia, a partir de los datos GEDI01_B)

• GEDI Level 2A: Producto de métricas de altura y elevación geolocalizadas (GEDI02_A).  El producto contiene 156 capas para cada una de las ocho haces. Tiene el propósito de proporcionar la interpretación de la forma de onda y productos extraídos de cada forma de onda recibida de GEDI01_B, incluidas las métricas de elevación del suelo, altura superior del dosel y altura relativa (RH). Para mayor información pueden consultar aquí, contando también con una guía de usuario en pdf.
  

Fuente: https://lpdaac.usgs.gov/resources/e-learning/getting-started-gedi-l2a-data-python/

• GEDI Level 2B: Producto de métricas de cobertura de dosel y perfil vertical (GEDI02_B), teniendo el propósito de extraer métricas biofísicas en cada forma de onda GEDI. Las métricas proporcionadas incluyen cobertura de dosel, índice de área de planta (PAI), densidad de volumen de área de planta (PAVD) y diversidad de altura de follaje (FHD). El producto contiene 96 capas para cada uno de los transectos terrestes de ocho haces. Para mayor información pueden ingresar aquí.

Fuente: https://lpdaac.usgs.gov/resources/e-learning/getting-started-gedi-l2b-data-python/

Descarga de datos lidar GEDI

Luego de conocer el tipo de datos en función a los productos disponibles, ahora vamos a seguir una secuencia de pasos para descargarlos, apoyándonos de las herramientas que nos proporcionan, entre ellas vamos a enfocarnos en el uso de Python.

Paso 1: Seleccionar nuestra área de interés, aquí es importante poder trabajarlo dentro de QGIS, primero podemos generar un ámbito de estudio en formato ESRI Shapefile y luego poder exportarlo a un formato geojson. Para nuestro caso se ha seleccionado una zona montañosa ubicada dentro del Departamento de Ucayali en Perú (ambito_pucallpa.geojson). Si vemos sus propiedades podemos contar con información de la extensión que ocupa en coordenadas.

Selección del área de interés

Paso 2: Usar el servicio GEDI Finder, el cual es un servicio web para localizar órbitas GEDI (archivos) que se cruzan con un cuadro delimitador de entrada. Para usarlo simplemente debemos adicionar los parámetros necesarios para encontrar los gránulos que contienen datos dentro de la región de interés del usuario, el mismo que ya definimos en el paso anterior. El esquema a seguir en caso queremos descargar un producto GEDI02_B para dicho área de interés es el que se muestra en la siguiente figura:

Modelo a seguir para emplear el servicio GEDI Finder

 Después de un tiempo en nuestro navegador web nos arroja el siguiente resultado:

Paso 3: Descargar directamente el archivo, para ello solo debemos hacer clic uno de los links que se generó, luego nos va a solicitar nuestras credenciales del EarthData, una vez ingresado nuestro usuario y contraseña la descarga empezará, es necesario guardarlo en una carpeta conocida. También se recomienda que la lista de links generados sean copiados a un editor de texto y guardarlo con la extensión *.csv.

Paso 4: Usar el script de Python GEDI_Subsetter, con la finalidad de convertir los productos de datos GEDI almacenados en formato HDF5, en archivos GeoJSON que se pueden cargar en SIG y software de percepción remota. Para ello debemos primero crear un "environment" en Anaconda, el mismo que debe contener algunas especificaciones en relación a la versión de Python y librerias requeridas. Simplemente nos vamos a un terminal de nuestro sistema para ingresar lo siguiente:

Creación de un ambiente con Conda denominado "gedi" junto a sus requisitos

Luego debemos activar el ambiente creado.

Activación del ambiente creado

Ahora podemos descargar el script GEDI_Subsetter.py dentro de nuestro ambiente de trabajo, una buena opción es clonarlo del repositorio en donde se ubica. Luego dentro de nuestro terminal ingresar la siguiente instrucción.

Ejecución del script de Python

Luego de un tiempo de procesamiento aparecerá dentro de la carpeta donde se ubica nuestro(s) productos GEDI, una subcarpeta denominada "output", en donde se guardará nuestro archivo en formato geojson.

Archivo en formato geojson generado

Paso 5: Visualización de nuestros datos en QGIS, puesto que ahora si lo podemos cargar para mostrar como podemos explorar los datos.

Vista luego de superponer un DEM y el archivo generado en formato geojson

En la imagen se puede apreciar las "huellas", las cuales están separadas a 60 metros y los transectos en distancias aproximadas de 600 metros.

Comprobación de las distancias de los transectos y de las huellas

Un paso necesario en abrir la tabla de atributos y en la columna que indica "l2b_quality_flag", eliminar los registros que tengan valores de cero.

Paso 6: Generar una visualización en 3D integrando una imagen y un archivo DEM, para ello vamos a necesitar tener instalado el plugin de QGIS denominado Qgis2threejs, el cual puede ser instalado desde el administrador de complementos. Para un mayor detalle de su uso puede ser consultado aquí.

Qgis2threejs dentro del administrador de complementos QGIS.

Ahora solo debemos activarlo y adicionamos a nuestro panel de capas una imagen de la zona que nos permitirá combinarla con nuestro DEM para generar una mejor visualización.

Capas listas para generar nuestra vista 3D

Debemos configurar nuestra capa "dem" para que se integre a nuestra imagen satelital (image6).

También debemos configurar nuestro capa del producto GEDI para tener una mejor visualización, como indicar que el tipo de objeto sea un cono, del mismo modo que el radio debe ser 25, considerando que la resolución espacial es de 25 metros y en la parte de altura, nuestra altura relativa (RH100) se divida por 100, con la finalidad de pasarlo de centímetros a metros. 

Si hacemos una ampliación veremos el detalle de las altitudes de los conos que simulan la geometría de la cobertura vegetal presente.

Muy bien, eso era lo que quería mostrarles, hay que tener en cuenta que el plugin te permite exportarlo en formatos de imagen y también en formato web, esta última opción es muy interesante.

Bueno, no quiero terminar sin antes agradecer los material de aprendizaje que brindar de manera libre los del programa ARSET de la NASA. Del mismo modo ellos han publicado un video que muestra gran parte del procedimiento seguido. 

Los datos empleados fueron del GEDI02_B:

Dubayah, R., Tang, H., Armston, J., Luthcke, S., Hofton, M., Blair, J. (2020). GEDI L2B Canopy Cover and Vertical Profile Metrics Data Global Footprint Level V001 [Data set]. NASA EOSDIS Land Processes DAAC. Accessed 2021-03-22 from https://doi.org/10.5067/GEDI/GEDI02_B.001

Referencias Consultadas:

1. GEDI Overview: https://lpdaac.usgs.gov/data/get-started-data/collection-overview/missions/gedi-overview/
2. GEDI-subsetter: https://git.earthdata.nasa.gov/projects/LPDUR/repos/gedi-subsetter/browse

Uso del Relief Visualization Toolbox (RVT) con Python y QGIS - Parte 2

 

Continuando con lo visto en la anterior entrada, ahora vamos a terminar de conocer los otros tipos de visualizaciones que podemos generar con la herramienta RVT. Se busca brindar el marco teórico sobre los que se sustentan las funciones y parámetors de la librería de Python que se emplean.

Factor de Vista del Cielo, Factor de Vista Anisotrópico y Apertura

Dentro del RVT, las tres visualizaciones mencionadas se calculan con la misma función, en el caso de la Apertura, la función calcula la Apertura positiva, mientras que en el caso de Apertura Negativa, solo debemos multiplar el DEM por -1 y con este DEM como entrada incluirla en la función. A continuación veremos algunos detalles que involucran a dichas visualizaciones.

Factor de Vista del Cielo - Sky-View Factor (svf)

Es una nueva técnica de visualización en relieve basada en iluminación difusa, en lugar de directa. Utiliza el factor de vista del cielo, un parámetro que corresponde a la porción de cielo visible limitada por el relieve. Utiliza el factor de vista del cielo como una técnica general de visualización del relieve para mostrar las características del relieve. En particular, mostramos que esta visualización es una herramienta muy útil en arqueología, ya que mejora el reconocimiento de características a pequeña escala desde un DEM de alta resolución. La SVF varía entre 0 y 1. Valores cercanos a 1 significan que casi todo el hemisferio es visible, como es el caso de las características expuestas (planos y picos), mientras que valores cercanos a 0 están presentes en sumideros profundos y partes bajas de valles profundos desde donde casi no se ve el cielo (1).

Fig. 1: El factor de vista del cielo se define por la parte del cielo visible (Ω) por encima de un cierto punto de observación visto desde una representación bidimensional (a). El algoritmo calcula el ángulo de elevación vertical del horizonte γi en n (aquí se presentan ocho) direcciones al radio especificado R (b). Fuente (1)

Factor de Vista del Cielo Anisotrópico - Anisotropic Sky-View Factor (asvf)

El factor de vista del cielo anisotrópico asume que el cielo es más brillante en algunas direcciones que en otras. Los pesos se basan en la función coseno de la mitad del ángulo (2). Hay tres parámetros que normalmente deben establecerse: el acimut del peso más alto, el exponente que define el gradiente del peso máximo al mínimo y el peso mínimo posible. El exponente y el peso mínimo definen el nivel de anisotropía: cuanto mayor es el exponente y menor el peso mínimo, más fuerte es el efecto. Esto se ha simplificado a los niveles de anisotropía bajo y alto donde los pesos mínimos posibles son 0.4 y 0.1, y los exponentes son 4 y 8 respectivamente (3).

Apertura - Openness (opns)

La apertura también es un sustituto de la iluminación difusa y se basa en una estimación de un ángulo medio de elevación del horizonte dentro de un radio de búsqueda definido. El valor medio de todos los ángulos cenitales da una apertura positiva, mientras que el valor medio del nadir da una apertura negativa (Fig. 2). La apertura positiva es similar al factor de vista del cielo, con una "sensación más aplanada", mientras que la apertura negativa brinda información adicional sobre las características convexas (4). Debido a que es independiente de la dirección y el sombreado y elimina la topografía general, es útil para el reconocimiento automático de características (3).

Fig. 2: Cálculo de la apertura positiva (α) y negativa (β) a lo largo de dos perfiles (a lo largo de la dirección Este (α 90, β90) y Oeste (α270, β270)) dentro de una distancia radial dada (r) en dos situaciones topográficas diferentes. Fuente (5). 

Uso de RVT

Para calcular las tres visualizaciones descritas debemos usar rvt.vis.sky_view_factor (). Los parámetros son: 

• dem (numpy.ndarray) - Ingrese el modelo de elevación digital como matriz numérica 2D.
• resolution (int) - Resolución de nuestro DEM
• compute_svf (bool)- Si es verdadero (True), calcula el factor de vista del cielo, 
• compute_asvf (bool)-  Si es verdadero (True), calcula svf anisotrópico, 
• compute_opns (bool)- Si es verdadero (True), calcula la apertura positiva), 
• svf_n_dir (int)- Número de direcciones, teniendo como rango valores de 4 a 360, pero el valor de 16 es óptimo para la mayoría de aplicaciones, 
• svf_r_max (int)-  Radio de búsqueda máximo en píxeles, recomendando valores pequeños como 5-10 si los elementos a distinguir son pequeños,
• svf_noise (int)- Nivel de eliminación de ruido (0-no eliminar, 1-bajo, 2-medio, 3-alto), 
• asvf_level(int)- Nivel de anisotropía (1-bajo, 2 -alto), 
• asvf_dir(int or float)- Dirección de anisotropía (0 es el Norte y 90 es el Este. 315 presenta la dirección habitual para el relieve sombreado), 
• ve_factor (int or float)- Factor de exageración vertical (Sin exageración = 1, para uso en relieve plano> 1),
• no_data (int or float)- Valor que representa no_data, todos los píxeles con este valor se cambian a np.nan.,
• fill_no_data (bool)- Si es Verdadero (True), se llena donde np.nan (no_data) con la media de píxeles circundantes (3x3),
• keep_original_no_data(bool)- . Si es Verdadero (True), cambia todos los píxeles de salida a np.nan donde dem tiene no_data.

La función genera un diccionario con claves: “svf” (si compute_svf es verdadero), “asvf” (si compute_asvf es verdadero), “opns” (si compute_opns es verdadero). Cada clave contiene una numpy array de visualización.

Uso del Relief Visualization Toolbox (RVT) con Python y QGIS - Parte 1

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:

1. Relief Visualization Toolbox in Python (https://rvt-py.readthedocs.io/en/latest/)
2. Understanding-azimuth-and-elevation (https://www.photopills.com/articles/understanding-azimuth-and-elevation)
3. Multi-Directional Hillshade Makes Your Maps Pop (https://www.esri.com/about/newsroom/arcuser/multi-directional-hillshade-makes-your-maps-pop/)
4. Función Sombreado (https://pro.arcgis.com/es/pro-app/latest/help/analysis/raster-functions/hillshade-function.htm)
5. Using  LIDAR-derived Local Relief Models (LRM) as a new tool for archaeological prospection (https://www.jstor.org/stable/j.ctt6wp79m.30)
6. RVT manual, version 2.2.1 (https://www.zrc-sazu.si/sites/default/files/rvt_2.2.1_0.pdf
7. 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)