Construyendo mapas de sombras hillshade con R

A raíz de la entrada 'Efecto hillshade o sombras de laderas' del blog de Dominic Royé, he probado una forma menos ambiciosa de obtener con R un sombreado hillshade de un mapa de elevaciones, pero con la premisa de no usar ninguna librería específica haciendo el método fácilmente trasladable a cualquier lenguaje con notación matricial como Python o Matlab. El código además está vectorizado para no usar bucles lo que hace su ejecución prácticamente instantánea.

Pese a su fuerte apariencia volumétrica un hillshade no es una infografía 3D, sino un cálculo en 2D de iluminación por zonas basado en el ángulo de incidencia de los rayos solares sobre el terreno.

Con un truco vectorial que ideé en la universidad estimamos la dirección normal al terreno en cada posición del mapa. Consiste en realizar el producto vectorial de dos vectores definidos por las elevaciones de las dos parejas de celdas adyacentes en los ejes X e Y. Elegante y sencillo, no requiere que las cuatro elevaciones participantes sean coplanarias.

Contrastando la dirección obtenida con la dirección de incidencia definida para la iluminación mediante un producto escalar, tenemos de forma directa el nivel de luminosidad reflejada en cada celda.

Se muestran las ecuaciones de los cálculos descritos incluyendo la forma de vectorizar la obtención del vector normal en todo el mapa (hacer clic para ampliar):

El sistema de coordenadas elegido hace coincidir la componente X con la dirección sur y la componente Y con la dirección este, siendo Z la elevación, y con este criterio implementaremos el vector iluminación en el código:

Aplicamos el cálculo a los datos DEM de la Sierra de Guadarrama en Madrid descargados desde el Centro Nacional de Información Geográfica del Ministerio de Transportes, que representados en escala de grises tienen el siguiente aspecto (hacer clic para visualizar en alta resolución):

El resultado del hillshade para una dirección de iluminación que correspondería al atardecer resulta del siguiente modo (hacer clic para visualizar en alta resolución):

Finalmente con Photoshop fusionamos la luminosidad del hillshade con el color del mapa de elevaciones coloreado con las curvas de la paleta 'topo.colors' que calculamos en 'Mapas de color 'viridis' en Photoshop con R' (hacer clic para ver en todo su esplendor):

Aplicando el mismo procedimiento a la impresionante orografía de Tenerife tenemos (hacer clic para alta resolución):

Y haciendo lo propio con la sierra de Alicante (clic para alta resolución):

Acabamos con un hillshade de Marte obtenido con datos de la NASA, marcando el punto de mayor y menor altitud de su superficie (clic para alta resolución, NOTA: son 32Mpx):

Hacemos zoom en el Olympus Mons, la mayor montaña conocida del Sistema Solar. Se trata de un volcán extinguido con una extensión cercana a Francia y una altitud respecto a su área circundante de más de 20km. Lo coloreamos con tonos marcianos y aplicamos un gradiente para simular el amanecer (procesado en Photoshop aquí, clic para mayor resolución):

~~~

En ejercicios posteriores que involucraron hillshades como 'Plegando papel con R', me di cuenta de que cuando necesitemos mezclar de algún modo dos mapas de elevación (ya sea promediándolos, por selección de zonas,...) para obtener un hillshade conjunto, es muy recomendable generar primero sus hillshades por separado y hacer a continuación la fusión de estos. El motivo es que fusionando mapas de elevación es muy probable que aparezcan pendientes bruscas o incluso discontinuidades en altitud que van a llenar de artefactos visibles el hillshade.

Otra reflexión es que pese a no ser un hillshade una verdadera infografía 3D como apuntábamos al principio, sí equivaldría al 100% a una representación tridimensional donde se aplicasen las siguientes restricciones:

• Dirección de observación cenital (de arriba abajo)
• Perspectiva axonométrica de la planta del terreno (sin fugas cónicas)
• Fuente(s) de iluminación según una dirección (fuente puntual en el infinito)
• Cálculo de iluminaciones de las caras visibles sin proyección de sombras (aunque esto lo paliamos en 'Proyección de sombras sobre un DEM con R' añadiendo al hillshade un cálculo preciso de sombras).

Cumpliendo todas esas condiciones cualquier render 3D proporcionaría una imagen indistinguible de un hillshade.

~~~

Repositorio con el código R: GitHub.