Dibuja tus mapas en R usando archivos de formato Shapefile

Este breve tutorial te mostrará como graficar mapas en R, extraer datos de archivos de formato shapefile, resaltar puntos en tu mapa y modificar su apariencia utilizando las bibliotecas rgdal y ggplot2.

Comenzando

Para que todas las bibliotecas funcionen correctamente, se recomienda instalar o actualizar a las últimas versiones de R y de su IDE RStudio con su correspondiente complemento RTools.

En este artículo, vamos a utilizar archivos de formato Shapefile. Si todavía no sabes de qué se trata, te explicamos a continuación.

Un archivo Shapefile contiene al menos:

• .shp - un archivo tipo shape, es la geometría misma.
• .shx - un archivo tipo index, tiene las posiciones indexadas del archivo .shp.
• .dbf - un archivo tipo attribute, tiene los atributos de cada forma en una columna, es de tipo dBase IV.

Adicionalmente, la carpeta donde se encuentran dichos archivos pueden contener otros archivos de formato .prj o .sbn, estos aportan más datos de la geometría o pueden ser usados en otros programas de sistemas de información geográfica.

Los datos importados con rgdal a partir de un archivo shapefile, son objetos que contienen polígonos vectoriales, con las coordenadas la latitud y la longitud en formato decimal. A partir de estos objetos, podemos extraer datos para graficarlos desde una tabla. La biblioteca broom nos ayuda en la extracción y agrupación de los datos para su correspondiente graficación.

Obteniendo nuestros archivos shapefile

Muchos de estos archivos Shapefile representan mapas de nuestros Estados, por lo que están disponibles de manera gratuita en la mayoría de los casos, en otros, son de paga o están más completos, actualizados y/o poseen datos específicos. Debajo mostramos algunos recursos web gratuitos que puedes usar.

data.humdata.org para shapefiles de varios países del mundo (algunos están desactualizados).

geo.gob.bo para shapefiles de Bolivia.

Instalación de pre-requisitos

Las siguientes bibliotecas de R son necesarias, para realizar nuestro ejemplo.

1install.packages("ggplot2")    # biblioteca para graficar.
2install.packages("rgdal")      # biblioteca para abrir archivo de formato shapefiles (datos geográficos)
3install.packages("broom")      # biblioteca usada para extraer datos del archivo importado de shapefiles
4install.packages("tidyverse")  # biblioteca usada para unir tablas
5install.packages("rio")        # biblioteca para importar datos de archivos csv
6install.packages("dplyr")     # biblioteca usada para agrupar por valores de columnas
7install.packages("extrafont") # biblioteca para importar los tipos de letra de windows.

(Recuerde que para ejecutar una línea de comando en el Editor de RStudio se usa Ctrl+Enter)

Preparación de los datos para graficar

Nuestra base de datos fue descargada de GeoBolivia, INE Bolivia y geodatos.net. Una vez descargados los datos, fueron depurados para el ejemplo.

Abrimos nuestras bibliotecas requeridas:

1library(ggplot2)
2library(rgdal)
3library(broom)
4library(tidyverse)
5library(rio)
6library(dplyr)
7library(extrafont)

Redireccionamos el directorio actual a nuestro directorio de trabajo e importamos nuestros archivos shapefiles:

1setwd("../mypath/") # redirecciona el directo actual a nuestro directorio de trabajo
2# importamos los datos geográficos a nuestra variable shapefile
3shapefile = readOGR(
4  dsn = ".",
5  layer = "departamentos_geo",
6  encoding = 'utf-8',
7  use_iconv = TRUE
8) 

• dsn: carpeta dentro del directorio actual, donde se encuentran ficheros shapefiles. Si se pone solo un punto "." hace referencia a que los ficheros se encuentran en la carpeta actual.

• encoding="utf-8", use_iconv=TRUE: indica a la función readOGR que debe importarse con la codificación utf-8. Como nuestro archivo contiene caracteres del español como ñ y vocales con tildes, es conveniente usar este comando.

Para observar el contenido de shapefile, use View(shapefile).

Extraemos datos en geotable y a continuación mostramos su cabecera:

1geotable = tidy(shapefile)
2head(geotable)

 1> head(geotable)
 2# A tibble: 6 x 7
 3   long   lat order hole  piece group id   
 4  <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>
 51 -65.8 -18.0     1 FALSE 1     0.1   0    
 62 -65.8 -18.0     2 FALSE 1     0.1   0    
 73 -65.8 -18.0     3 FALSE 1     0.1   0    
 84 -65.8 -18.0     4 FALSE 1     0.1   0    
 95 -65.8 -18.0     5 FALSE 1     0.1   0    
106 -65.8 -18.0     6 FALSE 1     0.1   0    

Nótese que geotable, no posee las etiquetas de los nombres de las regiones a graficar, esto lo arreglamos con:

1# añadiendo una columna id para poder juntar las columnas de nuestro geotable
2shapefile$id <- row.names(shapefile)
3
4# añadiendo a geotable los datos que faltan desde shapefile y lo juntamos por el id
5geotable <- left_join(geotable, shapefile@data, by = "id")

Para verificar que se han añadido los correspondientes nombres de regiones, use head(geotable)

Funciones auxiliares:

• shapefile$id <- : crea una nueva columna id en shapefile o la reemplaza.
• row.names(shapefile) : extrae los nombres de las filas del shapefile que por defecto es una numeración que va desde 0 y coincide con el id de nuestro geotable.
• shapefile@data accede a la tabla o dataframe data de nuestro shapefile.
• left_join(tabla1, tabla2, by = clave) : junta dos tablas por izquierda, es decir, añade valores que faltan de tabla2 a la tabla1 de acuerdo al código clave o columna común.

Ahora importamos nuestra tabla con los datos de población por departamento: (haciendo header = TRUE nos aseguramos que la primera fila se convierta en los nombres de las columnas)

1# importando la tabla de datos de población
2poblacion <- import("departamentospoblacion.csv", header = TRUE)

Nótese que cambiamos el nombre de nuestra columna en población DEPARTAMENTO por DEPARTAMEN, para que coincidan y luego juntarlas con left_join().

1# cambiando el nombre de la columna de DEPARTAMENTO POR DEPARTAMEN
2colnames(poblacion)[colnames(poblacion) == "DEPARTAMENTO"] <- "DEPARTAMEN"
3# juntamos en datos ambas tablas
4datos<-left_join(geotable, poblacion, by = "DEPARTAMEN")

Graficando el mapa con ggplot2

Ahora, graficamos con ggplot2:

1ggplot() +
2  geom_polygon(data = datos, aes(
3    x = long,
4    y = lat,
5    group = group,
6    fill = Poblacion2022
7  )) +
8  coord_equal() + labs(fill = "POBLACION")

• geom_polygon(data, aes(x, y, group, fill)): dibuja polígonos con la tabla data y muestra la estética con aes(). x, y son los valores en ejes de las abscisas y las ordenadas. Están agrupadas por el valor de la columna group y se colorea con base a los valores de la columna fill.
• coord_equal(): obliga a la gráfica a que la relación de aspecto entre coordenadas sea 1:1.
• lab(fill): pone el título a la leyenda con fill.

Mejorando la presentación de nuestro mapa

Entre las cosas que podemos hacer para mejorar la apariencia de nuestra gráfica están: poner un título, cambiar los colores, el fondo, el formato de nuestra leyenda. A continuación mostramos como hacerlo.

Extraemos los valores de la población para ponerlas como etiquetas dentro de nuestro mapa, para ello utilizamos las funciones group_by() y summarise() de la biblioteca dplyr.

1etiquetas_poblacion <- datos %>% group_by(DEPARTAMEN) %>%
2  summarise(
3    label_long = mean(range(long)),
4    label_lat = mean(range(lat)),
5    pob = mean(Poblacion2022)
6  )

• tabla0 %>% funcion0 %>% funcion1 ...: Esta notación indica que se debe tomar la tabla0 como argumento de la funcion0, luego los resultados de la funcion0 deben tomarse como argumentos de la función1 y así sucesivamente. El símbolo %>% es conocido como pipe operator, este nos ayuda a concatenar valores de entrada y salida de diferentes funciones.
• group_by(col) %>% summarise(col1=accion1, col2=accion2 ...): agrupa los datos en función del valor de columna col y con summarise() usa los datos agrupados para devolver nuevos valores: col1, col2 ..., que pueden estar en función de los valores agrupados.
• range(v): extrae los valores máximo y mínimo de un rango de datos v.
• mean(v): devuelve el valor promedio del vector v.

Si queremos incluir nuevas fuentes para el tipo de letra para nuestro mapa, empleamos los siguientes comandos: (Nota: nos pedirá confirmación para realizar la importación de fuentes y tardará unos minutos):

1extrafont::font_import("C:/Windows/Fonts")
2loadfonts(device = "win")
3fonts() # Nos muestra las fuentes disponibles.

Incluimos los datos de etiquetas_población en nuestra gráfica.

 1ggplot() +
 2  geom_polygon(data = datos, aes(
 3    x = long,
 4    y = lat,
 5    group = group,
 6    fill = Poblacion2022
 7  )) +
 8  coord_equal() + theme_void() +
 9  geom_text(
10    size = 4.5,
11    alpha = 0.9,
12    fontface = "bold",
13    data = etiquetas_poblacion,
14    mapping = aes(
15      x = label_long,
16      y = label_lat,
17      label = format(pob, big.mark = " "),
18      color = pob
19    )
20  ) +
21  labs(title = "Población estimada de Bolivia por departamentos para el año 2022",
22       fill = "Habitantes",
23       caption = "Datos:INE Boliva; GeoDatos: GeoBolivia") +
24  scale_colour_gradientn(colours = c("black", "black", "white", "white", "white"),
25                         guide = "none") +
26  scale_fill_continuous(
27    low = "#C4FFD1",
28    high = "#05693E",
29    guide = "colorbar",
30    labels = scales::label_number(big.mark = " ")
31  ) +
32  theme(
33    plot.title = element_text(
34      size = 14,
35      face = "bold",
36      family = "Helvetica",
37      hjust = 0.5
38    ),
39    legend.title = element_text(size = 12, family = "Rubik"),
40    plot.caption = element_text(family = "Helvetica")
41  )

• theme_void(): elimina el fondo y los ejes de nuestra gráfica.
• geom_text(size, alpha, fontface, data, mapping = aes(x, y, label), color): extrae los datos de data, para graficar el texto label en las coordenadas x e y, si se desea una diferenciación de colores se usa color. Con size, alpha y fontface, se establece el tamaño, la opacidad y la estética del texto respectivamente.
• format(v, big.mark): da el formato al valor v, indicando la separación de miles con big.mark (En nuestro ejemplo el separador de miles es solo el espacio " ").
• labs(title, fill, caption): con title, fill, caption pone el texto del título, la leyenda y el pie del gráfico respectivamente.
• scale_colour_gradientn(colours, guide): aplica una escala de colores a todos los valores asignados al argumento color. En nuestro ejemplo tenemos color dentro de la función geom_text(... aes(.. color = pob ...) ...), es decir, los valores de pob estarán coloreados según los valores de colours (la sintaxis color y colors pueden intercambiarse sin problema con colour y colours) y para que la guía de leyenda no se muestre usamos guide = "none".
• scale_fill_continuous(low, high, guide, labels): establece una escala de colores continua a los valores asignados a fill. En nuestro ejemplo fill está dentro de la función geom_polygon(... aes(... fill = Poblacion2022 ...) ...), es decir, los valores de la columna Poblacion2022 estarán afectados por esta función. Usamos low para el color correspondiente al valor más bajo y high para el color del valor más alto. Con guide = colorbar mostramos nuestra leyenda en forma de colorbar y conlabels modificamos la apariencia en la escala de nuestro colorbar.
• scales::label_number(big.mark=" "): usa la función label_number() de la biblioteca scales. Con esto modificamos la apariencia de los números de nuestra leyenda poniéndole un espacio " " como separador de miles. (Nótese que podemos usar scales::funcion() en lugar de library(scales) funcion()).
• theme(plot.title, legend.title, plot.caption): modifica la apariencia del título del gráfico, el título de la leyenda y el pie del gráfico respectivamente.
• element_text(size, face, family, hjust): extrae propiedades del texto para modificar el tamaño, la estética, el tipo y la posición en horizontal.

Añadiendo ubicaciones a nuestro mapa

Incluyendo las ubicaciones de ciudades capitales de departamentos a nuestro mapa:

 1ciudades = import("ciudades.csv") # importamos la localización de las ciudades capitales
 2
 3ggplot() +
 4  geom_polygon(
 5    data = datos,
 6    aes(
 7      x = long,
 8      y = lat,
 9      group = group,
10      fill = DEPARTAMEN
11    ),
12    color = "gray",
13    size = 0.5
14  ) +
15  geom_point(
16    alpha = 0.7,
17    data = ciudades,
18    mapping = aes(x = lat, y = long, colour = Ciudad),
19    size = 5
20  ) +
21  coord_equal() + labs(title = "Capitales de Departamento en Bolivia",
22                       color = "Ciudades Capitales",
23                       caption = "Fuente: geodatos.net, GeoBolivia") +
24  scale_fill_brewer(palette = 'PuBuGn', guide = "none") +
25  scale_color_manual(values = rainbow(9)) +
26  theme_void() +
27  theme(
28    plot.title = element_text(
29      size = 14,
30      face = "bold",
31      family = "Helvetica",
32      hjust = 0.5
33    ),
34    legend.title = element_text(size = 12, family = "Rubik"),
35    plot.caption = element_text(family = "Helvetica")
36  ) 

• geom_point(alpha, data, mapping = aes(x, y, colour), size): dibuja puntos dentro del gráfico, comparte similares argumentos con geom_poligon() y geom_text().
• scale_fill_brewer(palette, guide): similar a la función scale_fill_continuous(low, high, guide, labels) aplica una escala de colores tipo brewer a todos los objetos asignados a fill y con pallete seleccionamos el tipo de paleta de colores a aplicar.
• scale_color_manual(values): nos permite usar una escala de colores manual, values debe ser un vector que contenga los valores de los colores de la escala.
• raibow(9): devuelve un vector con 9 colores del arcoiris.

Asignando colores

Puedes asignar los colores simplemente usando su nombre en inglés. Para el blanco es white, para el rojo, red. También puedes utilizar el código hexadecimal, como #FF4500 para el rojo anaranjado; agruparlos en una escala de colores utilizando el comando c("red","#FF4500"...). Una página recomendable para seleccionar colores y obtener su código de color con un click es r-charts.com/colors/. También puedes emplear las funciones auxiliares que ofrece R, por ejemplo: scale_color/fill_brewer/viridis_. Estas proporcionan escalas predefinidas que podrían mejorar el impacto visual.

Guardando nuestro mapa

RStudio ofrece la posibilidad de exportar fácilmente desde su menú ubicado encima de la vista previa del gráfico "Export". Podemos optar por guardar nuestro mapa con mayor calidad o cierto formato y, para ello, podemos usar ggsave(), que nos permite exportar o guardar nuestro último gráfico ejecutado.

1ggsave(
2  filename = "grafica.png",
3  path = ".../mypath/",
4  scale = 1,
5  device = "png",
6  dpi = 320
7)

Guarda el mapa con el nombre filename en la ruta path, con la escala y formato de scale y device. Con dpi indicamos la cantidad de píxeles por pulgada, que es la calidad de nuestro archivo a exportar.

Referencias

• Repositorio del ejemplo
• Instituto Nacional de Estadística Bolivia
• GeoBolivia
• geodatos.net
• Robinlovelace
• Documentacion de R