Mapa de coropletas en ggplot2

Un mapa de coropletas es un tipo de mapa donde diferentes áreas geográficas se colorean en base a una variable asociada a esas áreas. Un mapa de coropletas proporciona una manera intuitiva de visualizar cómo una variable específica (como la densidad de población, los ingresos, etc.) varían a lo largo de diferentes áreas geográficas.

Usaremos los datos incluidos en el paquete giscoR package, que proporciona mapas en el formato sf y el conjunto de datos de ejemplo tgs00026, que incluye la renta disponible de los hogares privados a nivel NUTS2. En este ejemplo, usaremos el año 2016 como referencia.

Mapa base

En primer lugar, necesitaremos obtener los datos espaciales que contienen la información geográfica a ser usada en el gráfico. Usando giscoR, podemos seleccionar las geometrías NUTS con giscoR::gisco_get_nuts() y visualizar el objeto con geom_sf:

# install.packages("sf")
library(sf)
# install.packages("dplyr")
library(dplyr)
# install.packages("ggplot2")
library(ggplot2)
# install.packages("giscoR")
library(giscoR)

# Año
año_ref <- 2016

# Datos
nuts2 <- gisco_get_nuts(year = año_ref,
                        resolution = 20,
                        nuts_level = 2) %>%
         select(NUTS_ID, NAME_LATN)

# Mapa base
ggplot(nuts2) +
  geom_sf()

Límites y proyecciones

# Transformamos la forma
nuts2_3035 <- st_transform(nuts2, 3035)

ggplot(nuts2_3035) +
  geom_sf() +
  xlim(c(2200000, 7150000)) +
  ylim(c(1380000, 5500000))

Unir los datos y el mapa

# Filtro para seleccionar los datos de 2016
renta_disponible <- giscoR::tgs00026 %>%
  filter(time == año_ref) %>%
  select(-time)

nuts2_3035_datos <- nuts2_3035 %>%
  left_join(renta_disponible, by = c("NUTS_ID" = "geo"))

Mapa de coropletas básico

Como resultado del proceso anterior, tenemos disponible un objeto sf con los datos estadísticos a ser dibujados. Podemos proceder ahora a crear el mapa de coropletas básico con ggplot2, pasando los valores de interés (columna values) al argumento fill de aes:

ggplot(nuts2_3035_datos) +
  geom_sf(aes(fill = values)) +
  xlim(c(2200000, 7150000)) +
  ylim(c(1380000, 5500000))

ggplot(nuts2_3035_datos) +
  geom_sf(aes(fill = values), 
          color = "white",
          linetype = 1,
          lwd = 0.25) +
  xlim(c(2200000, 7150000)) +
  ylim(c(1380000, 5500000))

ggplot(nuts2_3035_datos) +
  geom_sf(aes(fill = values),
          color = NA) +
  xlim(c(2200000, 7150000)) +
  ylim(c(1380000, 5500000)) +
  scale_fill_viridis_c(option = "B")

Todavía podemos modificar más características del mapa para crear un gráfico más llamativo.

Mapa de coropletas avanzado

Vamos a modificar algunas características del mapa para crear uno más avanzado. Para ello usaremos puntos de corte y colores para crear un mapa categórico derivado de los datos continuos de la renta disponible. Además, vamos a dibujar países colindantes no incluidos en la muestra.

Empezaremos por obtener todos los países del mundo, que serán usados como el mapa base.

# Obtenemos todos los países y transformamos al mismo CRS
paises <- gisco_get_countries(year = año_ref,
                              resolution = 20) %>%
           st_transform(3035)

# Gráfico
choro_adv <- ggplot() +
  # Primera capa con todos los países
  geom_sf(data = paises, fill = "grey80", color = NA) +
  # Establece límites
  xlim(c(2200000, 7150000)) +
  ylim(c(1380000, 5500000))

choro_adv

Ahora, podemos empezar a preparar la capa de coropletas. Para ello eliminamos los valores con NA y creamos categorías discretas. También preparamos las etiquetas que serán incluidas en la leyenda del mapa.

# Limpiamos los NA del conjunto de datos original
nuts2_3035_datos_limpios <- nuts2_3035_datos %>%
  filter(!is.na(values))

# Cramos cortes y discretizamos los valores
br <- c(0, 10, 15, 20, 25, 30, 50) * 1000

nuts2_3035_datos_limpios$renta_disponible <- cut(nuts2_3035_datos_limpios$values,
                                                 breaks = br,
                                                 dig.lab = 5)

# Creamos etiquetas personalizadas - e.g. (0k-10k]
labs <- c(0, 10, 15, 20, 25, 30, 50)
labs_plot <- paste0("(", labs[1:6], "k-", labs[2:7], "k]")

Por último, elegimos la paleta Inferno incluida en la función hcl.colors (R >= 3.6) y procedemos a crear el mapa final:

# Paleta
pal <- hcl.colors(6, "Inferno", rev = TRUE, alpha = 0.7)

# Agregamos la capa del mapa base
choro_adv +
  # Agregamos la capa de las coropletas
  geom_sf(data = nuts2_3035_datos_limpios,
          aes(fill = renta_disponible), color = NA) +
  # Labs
  labs(title = "Renta disponible de los hogares privados",
       subtitle = "Regiones NUTS2 (2016)",
       caption = "Eurostat, © EuroGeographics para los límites administrativos",
       fill = "€") +
  # Paleta personalizada
  scale_fill_manual(values = pal,
                    drop = FALSE,
                    na.value = "grey80",
                    label = labs_plot,
                    # Legend
                    guide = guide_legend(direction = "horizontal",
                                         nrow = 1,
                                         label.position = "bottom")) +
  # Tema
  theme_void() +
  theme(plot.caption = element_text(size = 7, face = "italic"),
        legend.position = "bottom")