Mapa de coropletas en R

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. Este tipo de gráfico proporciona de manera intuitiva cómo una determinada variable (la población, los ingresos, …) varían entre las diferentes áreas.

En este tutorial utilizaremos el paquete giscoR para obtener el mapa de Italia en formato sf y usaremos el conjunto de datos de muestra de la misma librería llamado tgs00026, que proporciona la renta disponible de los hogares privados a nivel NUTS2. Para este tutorial utilizaremos el año 2013 como referencia.

Mapa base

En primer lugar, tendremos que 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 la función gisco_get_nuts y especificar el país que queremos dibujar dentro del argumento country. Para crear el mapa podemos aplicar la función st_geometry a nuestros datos y usar la función plot.

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

año_ref <- 2013

nuts2_IT <- gisco_get_nuts(
  year = año_ref,
  resolution = 20, 
  nuts_level = 2,
  country = "Italy") %>%
  select(NUTS_ID, NAME_LATN)

plot(st_geometry(nuts2_IT))

Proyección

# Proyección
nuts2_IT_32632 <- st_transform(nuts2_IT, 32632)

plot(st_geometry(nuts2_IT_32632))

Unir el mapa con los datos

Como los objetos geográficos en sf se comportan como data frames podemos hacr una unión con la función left_join para juntar los datos geográficos con los estadísticos.

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

nuts2_IT_32632_datos <- nuts2_IT_32632 %>%
  left_join(renta_disp, by = c("NUTS_ID" = "geo"))

Mapa de coropletas básico

Después de unir nuestros datos tendremos un objeto sf con los datos estadísticos a ser dibujados. Podemos proceder ahora a crear un mapa de coropletas básico dibujando la columna values del data frame.

plot(nuts2_IT_32632_datos[, "values"],
     breaks = "jenks",
     main = "Mapa de coropletas")

plot(nuts2_IT_32632_datos[, "values"],
     breaks = "jenks",
     nbreaks = 10,
     pal = hcl.colors(10),
     main = "Mapa de coropletas")

Mapa de coropletas avanzado

En esta sección vamos a personalizar algunas características del mapa anterior para crear un mapa más avanzado. Vamos a utilizar puntos de corte y colores personalizados para crear una variable categórica derivada de los datos de la renta disponible en Italia (datos continuos). Además, vamos a dibujar los países de alrededor no incluidos en la muestra.

En primer lugar, vamos a obtener todos los países del mundo a ser usados como mapa base, los transformamos al CRS de Europa y hacemos zoom sobre Italia.

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

# Gráfico
plot(st_geometry(paises),
     col = "grey80", border = NA,
     xlim = c(4600000, 4700000),
     ylim = c(1500000, 2600000))

En segundo lugar, comenzamos a preparar la capa con el mapa de coropeltas. Para ello limipamos los datos y creamos categorías discretas utilizando la función cut. También preparamos las etiquetas a ser incluidas en la leyenda del mapa final.

# Limpiamos los NA de los datos
nuts2_IT_32632_datos_clean <- nuts2_IT_32632_datos %>%
  filter(!is.na(values))

# Creamos puntos de corte y discretizamos valores
br <- c(10, 12, 14, 16, 18, 20) * 1000

nuts2_IT_32632_datos_clean$renta_disp <-
  cut(nuts2_IT_32632_datos_clean$values,
      breaks = br, dig.lab = 5)

# Creamos las etiquetas - e.g. (0k-10k]
labs <- c(10, 12, 14, 16, 18, 20)
labs_plot <- paste0("(", labs[1:5], "k-", labs[2:6], "k]")

Por último, utilizando la paleta Spectral de la función hcl.color procedemos a crear el mapa definitivo, uniendo el mapa base, agregando el mapa de coropletas de Italia y una leyenda:

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

# Mapa base
plot(st_geometry(paises),
  col = "grey80",
  border = NA,
  xlim = c(4600000, 4700000),
  ylim = c(1500000, 2600000),
  main = "Renta disponible de hogares privados")
mtext("Regiones NUTS2 de Italia (2013)", side = 3, adj = 0)

# Mapa de Italia
plot(st_transform(nuts2_IT_32632_datos_clean[, "renta_disp"],
                  crs = 3035), # Transformanos al CRS de Europa
     border = NA,
     pal = pal,
     categorical = TRUE,
     add = TRUE)

# Leyenda
legend(x = "topright",
       title = "€/hogar",
       cex = 0.9,
       bty = "n",
       inset = 0.02,
       y.intersp = 1.5,
       legend = labs_plot,
       fill = pal)