En este artículo vamos a practicar en R con mapas de elevaciones en formato matricial, leyendo con el paquete raster archivos GeoTIFF procedentes de General Bathymetric Chart of the Oceans (GEBCO). Este recurso permite descargar en formato GeoTIFF y con una alta resolución, el mapa de elevaciones de cualquier parte del planeta definida interactivamente, incluyendo tanto las zonas terrestres (altimetrías) como todos los fondos oceánicos (batimetrías).
Con esta información vamos a realizar un serie de animaciones que simulen la forma en que se configurarían los perfiles costeros de nuestro planeta si el nivel del mar comenzara a subir hasta cubrir los puntos más elevados de cada continente. En la parte final de cada vídeo haremos el ejercicio contrario, dibujando los bordes de costa que se tendrían si se retirasen las aguas hasta secarse por completo esa parte de la Tierra.
Quisiera hacer dos disclaimers sobre el ejercicio, para quisquillosos y curiosos:
- Las animaciones que vamos a hacer son ciencia ficción; se estima que no hay agua helada en todo el planeta como para hacer subir el nivel del mar ni 100m. Pero es que estamos simulando un Diluvio Universal de procedencia divina, y por ende inexplicable a ojos de la ciencia y la razón. Los vídeos no pretenden alertar del calentamiento global ni de ningún supuesto climático real, al menos más allá de los 100m de cubierta líquida.
- Cada píxel se dibuja del color del elemento (tierra o agua) que predomine en la superficie abarcada por dicho píxel en cada fotograma. Inevitable consecuencia del remuestreo hecho para adaptarnos a la resolución de salida, el píxel donde cae el Everest por ejemplo no toma como valor de elevación la de su pico (el máximo en la celda), sino un promedio del área circundante correspondiente, que lógicamente es menor. Por eso el píxel que contiene el Everest se "hunde" mucho antes de llegar a los 8.849m de agua, y lo mismo pasa con todas las cumbres célebres y las fosas más profundas.
Tras seleccionar la zona de interés en GEBCO y descargar el archivo GeoTIFF, desde R lo leemos y visualizamos con el paquete raster. La región descargada es superior a la zona de interés final porque se ha usado para varios ejercicios:
projectRaster()) a una proyección cónica conforme de Lambert, que comprime las zonas septentrionales expandiendo las más meridionales. Pasamos además de coordenadas en grados a km:
No existe la proyección perfecta, todas tienen ventajas e inconvenientes según la aplicación. Para el ejercicio que vamos hacer solo reproyectaremos Europa y Norteamérica, por tener superficies importantes alejadas del ecuador. Y dejaremos sin reproyectar Sudamérica, Asia y África. En este enlace hay una buena exposición de conceptos sobre proyecciones y GIS.
A este GeoTIFF de Europa le he dado varios usos, como por ejemplo este mapa que muestra en un mismo color rangos equivalentes de elevaciones y profundidades, para tomar noción de lo profundo que es el océano comparado con las altitudes emergidas (hacer clic para mayor resolución):
Pero para el Diluvio nos vamos a quedar solo con la parte sur del continente, con idea de no perder las Islas Canarias a la vez que logramos un formato apaisado. Así realizamos el siguiente recorte (función crop()):
Finalmente hacemos un reescalado del mapa (función resample()) a la resolución final deseada para el vídeo, que es mucho menor que la que nos proporciona el GeoTIFF de GEBCO. No nos preocupa demasiado preservar la relación de aspecto original sino adaptarnos a la resolución de salida buscada (en este caso Full HD: 1920 x 1080 píxeles):
Es en este remuestreo donde se "pierden" los mayores valores de altitudes y depresiones marinas, ya que implica un promediado de los valores altimétricos abarcados por cada píxel de salida. En la versión remuestreada la cota máxima es de 4.086m, cuando en el mapa previo al remuestreo teníamos hasta 4.613m y en los datos originales antes de reproyectar se tenían 4.655m, siendo ésa la cota en el centro de cada celda según el criterio de GEBCO (el mayor pico de la zona, el Mont Blanc, tiene 4.810m).
Las funciones usadas del paquete raster me han sorprendido por su extrema, y en mi opinión injustificada, lentitud. Con el remuestreo se da además la anomalía de que usar una interpolación bilineal resulte más rápida que un muestreo nearest neighbour, que debería ser casi instantáneo.
Consultado con su autor Robert Hijmans: "raster is no longer further developed. It has been replaced by terra and I think you will find that these methods are much faster in terra". Felicidades Guillermo por usar un paquete obsoleto.
En el GitHub de terra se ve uso intensivo de C++ así que tiene pinta de que las funciones críticas se han optimizado. Replicando con terra las operaciones del ejercicio tenemos mejoras dramáticas, especialmente en el remuestreo y visualización de los raster, aunque por desgracia no tanto en las reproyecciones:
~~~
Abandonamos la librería raster para hacer el resto de cálculos en R base, que van a ser todos muy rápidos al estar plenamente vectorizados.
Recogemos los datos de elevaciones en una matriz sobre la que obtenemos el histograma. En él ya podemos ver cómo los valores negativos (profundidades oceánicas) alcanzan valores mucho mayores a sus homónimas elevaciones sobre el nivel del mar:
Clasificamos los valores de este DEM matricial en tierra o mar según su valor de altitud esté por encima o por debajo de 0m (nivel del mar):
Al tener claramente definidas las categorías tierra/mar a nivel píxel, los contornos continentales se pueden obtener con precisión a partir de este mapa "sólido" con el simple filtrado paso alto detector de bordes que ya usamos (para otro fin) en 'Transferencia de calor por elementos finitos con R (I). Simulador'. Se hace en una línea de código matricial:
NOTA: los Grandes Lagos de Norteamérica, el Mar Caspio en Asia y los lagos Victoria, Tanganica y Malawi de África, tienen láminas de agua situadas en altitudes diferentes al nivel del mar. Para que sus contornos aparezcan dibujados tal y como los conocemos, he particularizado en ellos este cálculo al nivel actual de cada lago respecto al nivel del mar.
Ya solo nos queda obtener el hillshade del mapa de elevaciones, para lo cual usamos la rutina que definimos en 'Construyendo mapas de sombras hillshade con R'. Este gráfico, con el anterior contorno superpuesto, constituirá el fondo de todos los frames, alterando en él solo el color de cada píxel para diferenciar las zonas de tierra y mar en cada fotograma:
Con estos cálculos hechos, solo queda ejecutar un bucle para generar tantos fotogramas como necesitemos para lograr la duración impuesta por la banda sonora escogida para la animación. En cada fotograma se calculará un nivel de agua (que no seguirá una progresión lineal para que no se inunden rápidamente los continentes al principio), en base al cual se coloreará el mapa.
Para colorear cada píxel de cada fotograma, convertimos el hillshade desde su forma matricial en un array con 3 canales: rojo, verde y azul, modificando los correspondientes al color deseado para el píxel (en tierra añadimos rojo y quitamos azul, y en el mar quitamos rojo y además reducimos un poco el contraste para simular tierra sumergida).
Añadimos un rudimentario contador de metros que indicará el nivel del mar en cada fotograma, siendo 0m el nivel actual del mar. Como desconozco si hay alguna forma sencilla de "escribir" texto sobre una matriz numérica, lo he hecho a base de "palotes" dibujados con las primitivas gráficas matriciales creadas en 'Dibujando gráficos de mapa de bits con R'.
Los fotogramas generados desde R se convierten en vídeos por línea de comandos con FFmpeg, añadiéndoles el audio obtenido de YouTube.
No hay comentarios:
Publicar un comentario
Por claridad del blog, por favor trata de utilizar una sintaxis lo más correcta posible y no abusar del uso de emoticonos, mayúsculas y similares.