Proyección de sombras sobre un DEM con R

Dada la naturaleza de los mapas de elevación (DEM), donde en cada celda del terreno tenemos que preocuparnos de una única coordenada tridimensional, aunque con algunas restricciones podemos acercarnos de forma sencilla a un cálculo exacto de proyección de sombras, sin necesidad de recurrir a verdaderos algoritmos 3D complejos basados en trazado de rayos (los usados en software como Blender, rayshader,...).

El algoritmo que vamos a implementar no analiza sobre qué puntos del terreno proyecta sombras cada elevación sino que funciona al revés, para cada punto del mapa se evalúa si existe alguna elevación que bloquee por ocultación la luz incidente en él. De ser así consideramos que la posición se encuentra bajo la proyección de una sombra (hacer clic para ver a mayor tamaño):

Por simplicidad del código consideramos solamente valores de azimuth de la fuente iluminación desde los 4 puntos cardinales. Elegimos como DEM de pruebas la Sierra de Guadarrama, con datos de resolución 25m bajados del Centro Nacional de Información Geográfica:

Para un vector de iluminación definido por (X=0, Y=30, Z=5), que correspondería a un Sol cercano al amanecer de los equinoccios:

Se tiene el siguiente mapa de sombras (hacer clic para ver a su tamaño original):

Visto de esta manera es complicado interpretar la validez de las sombras, así que las superponemos sobre un hillshade del mapa. Haciendo clic puede verse un GIF animado donde se muestran en bucle los cálculos para diferentes ángulos de incidencia. Aunque el hillshade debería actualizarse con cada fuente de iluminación, se ha dejado fijo para resaltar exclusivamente el cambio en las sombras:

~~~

Se ha dejado el código abierto para extenderlo y que soporte cualquier azimuth de la fuente de iluminación, pero tal y como está podemos aplicar un pequeño workaround para lograr ángulos arbitrarios consistente en rotar el DEM previamente a realizar el cálculo de sombras. El resultado deberá rotarse de nuevo de vuelta en la misma cantidad.

Un ejemplo con rotación de 30º y posterior vector de iluminación (X=0, Y=40, Z=5). Rodeamos el DEM de negro lo que equivaldría al nivel del mar respecto a las elevaciones del mapa:

Superponiendo las sombras al hillshade del mapa rotado y luego deshaciendo la rotación para ambos, tenemos el resultado final para el ángulo buscado de 30º. Mantenemos la proyección de sombras más allá de los bordes del mapa para obtener cierto aspecto de diorama (hacer clic para ver en alta resolución):

~~~

Si al hillshade con sombras le sumamos una transformación de pseudo relieve escalonado tenemos una imagen con cualidades de auténtico render 3D: el hillshade proporciona el sombreado de laderas en función de la luz incidente, el escalonado del mapa de elevaciones genera el equivalente a una vista en volumen de tipo isométrico, y las sombras añaden realismo. El coloreado final y el bokeh se han aplicado en Photoshop (hacer clic para ver en alta resolución):

Aquí se nuestra el making of con la aportación individual de cada una de las tres funciones comentadas.

~~~

Para terminar añadimos, con mínimo esfuerzo en código, estilos de sombra adicionales al modelo base de sombra "dura". Cuantifican según diferentes variables el grado de ocultación que sufre cada punto sombreado, produciendo sombras difusas más realistas (hacer clic para ver en alta resolución):

• 'hard': las sombras booleanas estándar que hemos usado hasta ahora.
• 'width': sombreado más intenso cuanto mayor sea el grosor de la zona que proyecta la sombra.
• 'max': sombreado más intenso cuanto mayor sea la diferencia de cota máxima entre el rayo de luz y la elevación que lo obstruye.
• 'mixed': combinación de los dos anteriores, da las sombras más suaves.

Se muestra aquí el resultado de cada tipo de sombreado aplicado a un relieve sintético (gorrito de Gandalf). A la izquierda sombras sobre el hillshade, a la derecha sobre el mapa de relieve escalonado (hacer clic para ver a mayor resolución):

~~~

Más allá del grafismo, la detección de sombras sobre el terreno tiene múltiples aplicaciones, desde poder realizar cálculos de insolación (cultivos, plantas fotovoltaicas,...), como también para conocer de antemano los ángulos admisibles en un barrido LIDAR para no perder por ocultación física datos en la captura.

En cuanto a la implementación he vectorizado parcialmente el código con Vectorize(), pero parece que esta función solo vectoriza a efectos formales de sintaxis del código sin optimizar realmente la ejecución: va más lento que recorriendo todo el mapa digital con dos bucles anidados, que es lo que he acabado haciendo. En definitiva éste ha sido un caso de vectorización fallida.

Repositorio con el código R y archivos DEM: GitHub.