Cartografía firefly

Cartografía firefly

Los mapas firefly (lampíridos ingl.) son promocionados y descritos por John Nelson quien publicó un post en 2016 sobre sus características. No obstante, este tipo de mapas están vinculados a ArcGIS, lo que me ha llevado a intentar recrearlos en R. La reciente extensión de ggplot2 con el paquete ggshadow nos facilitará la creación de este estilo cartográfico. Se caracteriza por tres elementos 1) un mapa base oscuro y desaturado (p.j. imágenes satelitales) 2) una viñeta y área resaltada enmascarada y 3) una única capa temática brillante. Lo esencial son los colores y el brillo que se logra con colores fríos, habitualmente colores neón. John Nelson explica más detalles en este post.

¿Para qué sirve el estilo firefly? En palabras de John Nelson: “the map style that captures our attention and dutifully honors the First Law of Geography”. John hace referencia a lo dicho por Waldo Tobler “everything is related to everything else, but near things are more related than distant things” (Tobler 1970).

En este post visualizaremos todos los terremotos registrados en el suroeste de Europa con una magnitud mayor de 3.

Paquetes

Usaremos los siguientes paquetes:

Paquete Descripción
tidyverse Conjunto de paquetes (visualización y manipulación de datos): ggplot2, dplyr, purrr,etc.
plotwidgets Contiene funciones para la conversión de colores (RGB, HSL)
terra Importar, exportar y manipular raster (paquete sucesor de raster)
raster Importar, exportar y manipular raster
sf Simple Feature: importar, exportar y manipular datos vectoriales
ggshadow Extensión de ggplot2 para geometrías con sombreado
ggspatial Extensión de ggplot2 para objetos espaciales
ggnewscale Extensión de ggplot2 para crear multiples scalas
janitor Funciones sencillas para examinar y limpiar datos
rnaturalearth Mapas vectoriales del mundo ‘Natural Earth’
# instalamos los paquetes si hace falta

if(!require("tidyverse")) install.packages("tidyverse")
if(!require("sf")) install.packages("sf")
if(!require("terra")) install.packages("terra")
if(!require("raster")) install.packages("raster")
if(!require("plotwidgets")) install.packages("plotwidgets")
if(!require("ggshadow")) install.packages("ggshadow")
if(!require("ggspatial")) install.packages("ggspatial")
if(!require("ggnewscale")) install.packages("ggnewscale")
if(!require("janitor")) install.packages("janitor")
if(!require("rnaturalearth")) install.packages("rnaturalearth")

# paquetes

library(raster)
library(terra)
library(sf)
library(tidyverse)
library(plotwidgets)
library(ggshadow)
library(ggspatial)
library(ggnewscale)
library(janitor)
library(rnaturalearth)

Preparación

Datos

Primero descargamos todos los datos necesarios. Para el mapa de base usaremos la Blue Marble vía el acceso a worldview.earthdata.nasa.gov donde me he descargado una selección del área de interés en formato geoTiff con una resolución de 1 km. Es importante ajustar la resolución al detalle necesario del mapa.

Importar

Lo primero que hacemos es importar el raster RGB Blue Marble y los datos de los terremotos. Para importar el raster hago uso del nuevo paquete terra que es el sucesor del paquete raster. Podéis encontrar una comparación reciente aquí. No todos los paquetes son compatibles todavía con la nueva clase de SpatRaster, por eso, nos hace falta también el paquete raster.

# terremotos

terremotos <- read.csv2("catalogoComunSV_1621713848556.csv")
str(terremotos)
## 'data.frame':    149724 obs. of  10 variables:
##  $ Evento       : chr  "          33" "          34" "          35" "          36" ...
##  $ Fecha        : chr  "  02/03/1373" "  03/03/1373" "  08/03/1373" "  19/03/1373" ...
##  $ Hora         : chr  "    00:00:00" "    00:00:00" "    00:00:00" "    00:00:00" ...
##  $ Latitud      : chr  "     42.5000" "     42.5000" "     42.5000" "     42.5000" ...
##  $ Longitud     : chr  "      0.7500" "      0.7500" "      0.7500" "      0.7500" ...
##  $ Prof...Km.   : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ Inten.       : chr  "     VIII-IX" "            " "            " "            " ...
##  $ Mag.         : chr  "            " "            " "            " "            " ...
##  $ Tipo.Mag.    : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ Localización: chr  "Ribagorça.L" "Ribagorça.L" "Ribagorça.L" "Ribagorça.L" ...
# Blue Marble RGB raster

bm <- rast("snapshot-2017-11-30T00_00_00Z.tiff")
bm # contiene tres capas (red, green, blue)
## class       : SpatRaster 
## dimensions  : 7156, 7156, 3  (nrow, ncol, nlyr)
## resolution  : 0.008789272, 0.008789272  (x, y)
## extent      : -33.49823, 29.39781, 15.77547, 78.67151  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : snapshot-2017-11-30T00_00_00Z.tiff 
## colors RGB  : 1, 2, 3 
## names       : snapshot-2~0_00_00Z_1, snapshot-2~0_00_00Z_2, snapshot-2~0_00_00Z_3
# plot

plotRGB(bm)

# límites países

limits <- ne_countries(scale = 50, returnclass = "sf")

Terremotos

En este paso limpiamos los datos importados de los terremotos. 1) Convertimos en númerico longitud, latitud y magnitud usando la función parse_number() y limpiampos los nombres de las columnas con la función clean_names(), 2) Creamos un objeto espacial sf y lo proyectamos usando el EPSG:3035 correspondiendo a ETRS89-extended/LAEA Europe.

# limpiamos los datos y creamos un objeto sf

terremotos <-  terremotos %>% clean_names() %>%
               mutate(across(c(mag, latitud, longitud),                                                                                                  parse_number)) %>%
               st_as_sf(coords = c("longitud", "latitud"), 
                       crs = 4326) %>% 
               st_transform(3035) # proyectamos a Laea 

Mapa de fondo Blue Marble

Recortamos el mapa de fondo a una extensión menor, pero todavía no limitamos el área final.

# recortamos al área deseada

bm <- crop(bm, extent(-20, 10, 30, 50))

Para obtener una versión desaturada del raster RGB de Blue Marble, debemos aplicar una función creada con este fin. En ésta usamos la función rgb2hsl() del paquete plotwidgets, que nos ayuda en convertir RGB a HSL y viceversa. El modelo HSL se define con Hue (tono), Saturation (saturación), Lightness (luminosidad). Los últimos dos parametros se expresan en ratio o porcentaje. El tono está definido en una rueda de color de 0 a 360º. 0 es rojo, 120 es verde, 240 es azul. Para cambiar la saturación únicamente debemos bajar el valor de S.

# función para cambiar la saturación desde RGB

saturation <- function(rgb, s = .5){
  
  hsl <- rgb2hsl(as.matrix(rgb))
  hsl[2, ] <- s
  
  rgb_new <- as.vector(t(hsl2rgb(hsl)))
  
  return(rgb_new)
  
}

Empleamos nuestra función saturation() usando otra función app() que la aplica a cada píxel con las tres capas de RGB. Añadimos el argumento s, el que define el nivel de saturación deseada. Este paso puede tardar varios minutos. Después proyectamos nuestra imagen RGB.

# aplicamos la función para desaturar con un 5%

bm_desat <- app(bm, saturation, s = .05)

# plot nuevo imagen RGB

plotRGB(bm_desat)

# proyectamos 

bm_desat <- terra::project(bm_desat, "epsg:3035")

Construcción del mapa firefly

Límites y gradícula

Antes de empezar a construir el mapa, creamos una gradícula y estabelecemos la extensión final del mapa.

# definimos los límites finales del mapa

bx <- tibble(x = c(-13, 6.7), y = c(31, 47)) %>% 
       st_as_sf(coords = c("x", "y"), crs = 4326) %>%
        st_transform(3035) %>% 
         st_bbox()

# crearmos una gradícula del mapa

grid <- st_graticule(terremotos) 

Mapa con fondo de imagen

La función layer_spatial() de ggspatial nos permite añadir un raster RGB sin grandes problemas, no obstante, todavía no apoya la nueva clase SpatRaster. Por eso, debemos convertirlo en la clase stack con la función stack(). También es posible usar en lugar de geom_sf(), la función layer_spatial() para objetos vectoriales de clase sf o sp.

ggplot() +
  layer_spatial(data = stack(bm_desat)) + # mapa de fondo Blue Marble
  geom_sf(data = limits, fill = NA, size = .3, colour = "white") + # límites países
  coord_sf(xlim = bx[c(1, 3)], 
           ylim = bx[c(2, 4)], 
           crs = 3035,
           expand = FALSE) +
  theme_void()

Mapa con fondo y los terremotos

Para crear el efecto del brillo en los mapas firefly, hacemos uso de la función geom_glowpoint() del paquete ggshadow. También existe la misma función para líneas. Dado que nuestros datos son espaciales de clase sf y la geometría no apoya directamente el uso de sf, debemos indicar como argumento stats = "sf_coordinates" y dentro de aes() indicar geometry = geometry. Mapearemos el tamaño de los puntos en función de la magnitud. Además, filtramos aquellos terremotos con una magnitud mayor del 3.

Dentro de la función geom_glowpoint(), 1) definimos el color deseado para el punto y el brillo, 2) el grado de transparencia con alpha bien para el punto o bien para el brillo. Por último, en la función scale_size() estabelecemos el rango (mínimo, máximo) del tamaño que tendrán los puntos.

ggplot() +
  layer_spatial(data = stack(bm_desat)) +
  geom_sf(data = limits, fill = NA, size = .3, colour = "white") +
  geom_sf(data = grid, colour = "white", size = .1, alpha = .5) +
  geom_glowpoint(data = filter(terremotos, mag > 3),
                 aes(geometry = geometry, size = mag), 
                   alpha = .8,
                   color = "#6bb857",
                   shadowcolour = "#6bb857",
                   shadowalpha = .1,
                   stat = "sf_coordinates",
                   show.legend = FALSE) +
  scale_size(range = c(.1, 1.5)) +
  coord_sf(xlim = bx[c(1, 3)], 
           ylim = bx[c(2, 4)], 
           crs = 3035,
           expand = FALSE) +
  theme_void()

Mapa final

El brillo de los mapas firefly se caracteriza por tener un tono blanco o un tono más claro en el centro de los puntos. Para lograrlo, debemos duplicar la capa anterior creada, cambiando únicamente el color y hacer los puntos con su brillo más pequeños.

Por defecto, ggplot2 no permite emplear diferentes escalas para la misma caracteristica (tamaño, color, etc) de distintas capas. Pero el paquete ggnewscale nos da la posibilidad de incorporar múltiples escalas de una característica de distintas capas. Lo único importante para lograrlo es el orden en el que se añade cada capa y escala. Primero debemos poner la geometría y despúes su escala correspondiente. Indicamos con new_scale('size') que la siguiente capa y escala es una nueva independiente de la anterior. Si usaramos colour o fill sería con new_scale_*().

ggplot() +
  layer_spatial(data = stack(bm_desat)) +
  geom_sf(data = limits, fill = NA, size = .3, colour = "white") +
  geom_sf(data = grid, colour = "white", size = .1, alpha = .5) +
  geom_glowpoint(data = filter(terremotos, mag > 3),
                   aes(geometry = geometry, size = mag), 
                   alpha = .8,
                   color = "#6bb857",
                   shadowcolour = "#6bb857",
                   shadowalpha = .1,
                   stat = "sf_coordinates",
                   show.legend = FALSE) +
  scale_size(range = c(.1, 1.5)) +
  new_scale("size") +
  geom_glowpoint(data = filter(terremotos, mag > 3),
                   aes(geometry = geometry, size = mag), 
                   alpha = .6,
                   shadowalpha = .05,
                   color = "#ffffff",
                   stat = "sf_coordinates",
                   show.legend = FALSE) +
  scale_size(range = c(.01, .7)) +
  labs(title = "TERREMOTOS") +
  coord_sf(xlim = bx[c(1, 3)], ylim = bx[c(2, 4)], crs = 3035,
           expand = FALSE) +
  theme_void() +
  theme(plot.title = element_text(size = 50, vjust = -5, colour = "white", hjust = .95))

ggsave("firefly_map.png", width = 15, height = 15, units = "in", dpi = 300)

Buy Me A Coffee

Dr. Dominic Royé
Dr. Dominic Royé
Investigador postdoctoral senior

Relacionado