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.

• Selección de Blue Marble via worldview.earthdata.nasa.gov (~66 MB) descarga
• Registros de terremotos históricos en el suroeste de Europa del IGN descarga

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 colorize(), 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.

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

bm_desat <- colorize(bm, to = "hsl")
bm_desat[[2]] <- .05 # ratio de saturación
set.RGB(bm_desat, 1:3, "hsl")
bm_desat <- colorize(bm_desat, to = "rgb")

# plot nuevo imagen RGB

plotRGB(bm_desat)

# proyectamos 

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

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)