Más sobre ciencia de datos: cienciadedatos.net

Introducción


En este documento se muestra cómo crear mapas en R empleando ggplot2. A diferencia de otros aproximaciones, en las que el mapa es una imagen sobre la que se superponen elementos, en este caso todo el mapa se crea empleando las geometrías de ggplot2. Además, una vez creado el gráfico, puede convertirse en interactivo con el paquete plotly.

A modo de ejemplo, se representan todos los municipios de España (la península y Baleares) y se colorea en función del número de habitantes del censo 2019.

Librerías


Las librerías utilizadas a lo largo del documento son:

library(tidyverse)
library(stringr)
library(rgdal)
library(rgeos)
library(plotly)



Mapa


Para crear un mapa es necesario disponer de las coordenadas que definen los límites de cada elemento, en este caso los municipios. Uno de los estándares para almacenar y compartir elementos geográficos y los atributos asociados a ellos es el formato Shapefile (.shp).

Se pueden descargar los ficheros .shp con todos los municipios españoles en las páginas del © Centro Nacional de Información Geográfica (CNIG) y del © Instituto Geográfico Nacional de España.

Datos adicionales


El código de identificación de las provincias y comunidades autónomas de España puede encontrarse en la página del INE. Instituto Nacional de Estadística.

El número de habitantes a nivel de comunidad autónoma, provincia y municipio, desde 1999 a 2020 puede descargarse de la página del padrón oficial del INE. Es recomendable descargarse todos los años en un único archivo .zip link.

Lectura de datos .shp


Los ficheros .shp pueden leerse en R con la función readOGR() del paquete rgdal.

# Lectura del fichero .shp
mapa <- rgdal::readOGR(
          paste0("./datos/lineas_limite/recintos_municipales_inspire_peninbal_etrs89/",
                 "recintos_municipales_inspire_peninbal_etrs89.shp")
        )

El objeto devuelto por readOGR() es una clase S4 por lo que sus elementos (llamados slots) son accesibles mediante el operador @. En su interior contiene, entre otros elementos, las coordenadas de cada polígono @polygons e información adicional sobre el mapa @data.

# Nombre de los slots del objeto mapa
slotNames(mapa)
## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"



Conversión de formato


Una vez que el fichero se han leído, para poder emplear las magníficas funcionalidades gráficas de ggplot2, se tiene que convertir a formato dataframe o tibble. Esta conversión puede hacerse con la función fortify() del paquete fortify. El argumento región es el identificador único empleado para asociar cada punto a un polígono. Con la función names() se pueden ver cuales son las opciones disponibles.

names(mapa)
## [1] "INSPIREID"  "COUNTRY"    "NATLEV"     "NATLEVNAME" "NATCODE"   
## [6] "NAMEUNIT"   "CODNUT1"    "CODNUT2"    "CODNUT3"
mapa_df <- fortify(model = mapa, region = "NATCODE")
mapa_df %>% head()



Metadatos


Además de las coordenadas, dentro del objeto SpatialPolygonsDataFrame hay información adicional sobre el mapa.

info_municipios <- mapa@data
info_municipios %>% head()

La variable NAMEUNIT es el nombre de la unidad, en este caso el nombre del municipio. La variable NATCODE es un código de identificación que sigue la siguiente estructura:

  • NATCODE: 34074000000
  • 34: País
  • 07: Comunidad Autónoma (Castilla y León)
  • 40: Provincia (Segovia)
  • 00000: Código Municipio

Se crean las columnas pais, c_autonoma, provincia y municipio a partir del NATCODE.

info_municipios <- info_municipios %>%
                    mutate(
                      pais       = str_sub(string = NATCODE, start = 1, end = 2),
                      c_autonoma = str_sub(string = NATCODE, start = 3, end = 4),
                      provincia  = str_sub(string = NATCODE, start = 5, end = 6),
                      municipio  = str_sub(string = NATCODE, start = 7, end = -1)
                    ) %>%
                    rename(nombre_municipio = NAMEUNIT)

# Se seleccionan las columnas de interés
info_municipios <- info_municipios %>%
                   select(
                     NATCODE, nombre_municipio, c_autonoma, provincia, municipio
                   )

Para poder unir la información de los municipios con el dataframe creado mediante fortify, es necesario un id común. En este caso, como se ha indicado region = "NATCODE" esta variable es el id común. Si no se hubiese especificado región, el id se corresponde con el orden de las filas (cuando se crea el dataframe con fortify(), el id se inicia en 0).

# Se añade la información de los municipios
mapa_df <- mapa_df %>%
           left_join(info_municipios, by = c("id" = "NATCODE"))



Gráfico


Cuantos más puntos (coordenadas) se emplean para crear el mapa, mayor es su resolución, pero también mayor el tiempo necesario para crear. Es este caso se reduce la resolución empleando únicamente 1 de cada 5 puntos.

# Se eliminan puntos (se reduce la resolución)
mapa_df <- mapa_df %>%
           slice(seq(1, nrow(mapa_df), 5))
mapa_df %>%
ggplot(aes(x = long, y = lat, group = group)) +
  geom_polygon(color = "gray20", fill = "white") + 
  coord_map("mercator") +
  theme_bw() +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )



Color


Al tratarse de un gráfico ggplot, se puede colorear cada elemento en función de cualquiera de las columnas del dataframe que contiene los datos del mapa.

Color por comunidad autónoma

mapa_df %>%
ggplot(aes(x = long, y = lat, group = group, fill = c_autonoma)) +
  geom_polygon(color = "black") +
  coord_map("mercator") +
  labs(title = "Municipios de España (peninsula y Baleares)",
       subtitle = "Color por comunidad autónoma") +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

Resaltar zonas

ggplot() +
  geom_polygon(
    data = mapa_df,
    aes(x = long, y = lat, group = group),
    fill = "gray90",
    color = "gray20") +
  geom_polygon(
    data = mapa_df %>% filter(c_autonoma %in% c("10", "12")),
    aes(x = long, y = lat, group = group),
    fill = "firebrick",
    alpha = 0.7,
    color = "black") +
  coord_map("mercator") +
  labs(title = "Municipios de España (peninsula y Baleares)",
       subtitle = "Coloreada la Comunidad Valenciana y Galicia") +
  theme_bw() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
)

Color por número de habitantes

Se colorea cada municipio en función del número de habitantes.

# Datos del censo municipal de 2019
pob_municip_19 <- readxl::read_excel(
                      path = "./datos/pobmun19.xlsx",
                      skip = 1
                    )

head(pob_municip_19)

Para tener una idea de la magnitud de las cifras, se muestran los 5 municipios con mayor y menor población.

pob_municip_19 %>% arrange(desc(POB19)) %>% head(5)
pob_municip_19 %>% arrange(POB19) %>% head(5)

Para poder colorear el mapa en función de la población, hay que añadir la columna POB19 a mapa_df. Afortunadamente, esta unión puede hacerse empleando el código NATCODE en lugar del nombre de los municipios (evitando así problemas de acentos y otros caracteres que dificultan la coincidencia).

La columna municipio del mapa_df se corresponden con la unión de las columnas CPRO y CMUN de pob_municip_2019.

pob_municip_19 <- pob_municip_19 %>% 
                   mutate(
                      municipio = stringr::str_c(CPRO, CMUN)
                   )
mapa_df <- mapa_df %>% 
           left_join(
              y  = pob_municip_19 %>% select(municipio, POB19),
              by = "municipio"
           )
mapa_df %>%
ggplot(aes(x = long, y = lat, group = group, fill = log10(POB19))) +
  geom_polygon(color = "black") +
  coord_map("mercator") +
  scale_fill_viridis_c() +
  labs(title = "Municipios de España (peninsula y Baleares)",
       subtitle = "Color por población, Escala logarítmica") +
  theme_bw() +
  theme(
    legend.position = "bottom",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
)



Regiones del mapa


Para representar únicamente una región, se puede filtrar el dataframe por cualquiera de las columnas pero, para que tenga sentido el mapa, no de deben de romper las unidades que forman los polígonos, en este caso los municipios.

# Se selecciona la comunidad autónoma
mapa_df_c_valenciana <- mapa_df %>% filter(c_autonoma == "10")
mapa_df_c_valenciana %>%
ggplot(aes(x = long, y = lat, group = group, fill = log10(POB19))) +
  geom_polygon(color = "black") +
  coord_map("mercator") +
  scale_fill_viridis_c() +
  labs(title = "Municipios de la Comunidad Valenciana",
       subtitle = "Color por población, Escala logarítmica") +
  theme_bw() +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

# Se selecciona Valencia
mapa_df_valencia <- mapa_df %>% filter(provincia == "46")
mapa_df_valencia %>%
ggplot(aes(x = long, y = lat, group = group, fill = log10(POB19))) +
  geom_polygon(color = "black") +
  coord_map("mercator") +
  scale_fill_viridis_c() +
  labs(title = "Municipios de Valencia",
       subtitle = "Color por población, Escala logarítmica") +
  theme_bw() +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )



Etiquetar el mapa



Con la función coordinates() identifica el centro (aproximado) de cada polígono, lo que perite saber dónde añadir el texto.

nombres_municipios <- mapa@data$NAMEUNIT
centro_municipios  <- coordinates(mapa)
colnames(centro_municipios)  <- c("long", "lat")
centro_municipios  <-  as.data.frame(centro_municipios)
centro_municipios$nombre_municipio <- nombres_municipios
ggplot() +
  geom_polygon(
    data = mapa_df %>% filter(provincia == "41"),
    aes(x = long, y = lat, group = group),
    fill = "gray90",
    color = "black") + 
  geom_polygon(
    data = mapa_df %>% filter(nombre_municipio == "Écija"),
    aes(x = long, y = lat, group = group),
    fill = "firebrick",
    alpha = 0.5) + 
  geom_text(
    data = centro_municipios %>% filter(nombre_municipio == "Écija"),
    aes(x = long, y = lat, label = nombre_municipio),
    inherit.aes = FALSE,
    size = 8
  ) +
  coord_map("mercator") +
  theme_bw() +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks =  element_blank(),
    axis.title = element_blank(),
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )



Mapa interactivo


p <- mapa_df_c_valenciana %>%
    ggplot(aes(x = long, y = lat, group = group, fill = log10(POB19),
               texto = nombre_municipio)) +
      geom_polygon(color = "black") +
      coord_map("mercator") +
      scale_fill_viridis_c() +
      labs(title = "Municipios de Valencia",
           subtitle = "Color por población, Escala logarítmica") +
      theme_bw() +
      theme(
        legend.position = "none",
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks =  element_blank(),
        axis.title = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()
      )

ggplotly(p, width = 450, height = 700, tooltip = c("texto"))



¿Cómo citar este documento?

Mapa de España con ggplot2 y R por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/58_mapas_con_r.html



Creative Commons Licence
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.