Orlando Sabogal

Now researching in cycling and social capital

Hagamos un mapa: librería tmap

Orlando Sabogal / 2019-07-23

Como probablemente a muchas otras personas les ha pasado, la primera vez que vi La Marcha de Napoleón de 1812-1813 elabarado por Charles Joseph Minard en 1869 quedé impresionado. El mapa (ver siguiente Figura) narra la historia de la fallida campaña Napoleónica a Rusia y su regreso a Francia. La invasión comenzó el 12 de Marzo y durante mucho tiempo el ejército Francés avanzó sin que las fuerzas Rusas pudieran brindar mucha resistencia. Entre el 2 y 6 de Septiembre los generales Rusos ordenaron evacuar Moscú y quemar la ciudad, por lo que Napoleón encontró una ciudad sin las condiciones para albergar a sus soldados en invierno. Esto, sumado a que Rusia estaba reorganizando su ejército llevó a que el 19 de Octubre los franceses iniciaran la retirada. Pero la retirada tampoco sería fácil, el invierno y los Rusos los persiguieron sin piedad.

Lo que hace tan interesante el mapa de Minard es que logra encapsular mucha información en una sola imagen: cantidad de soldados, dirección, ubicación (longitud y latitud), división de las tropas, reagrupramiento de las tropas y temperatura. Edward Tufte en su libro The Visual Display of Quantitative Information lo resalta como el mejor gráfico estadístico jamás realizado.

Campaña Napoleónica contra Rusia en 1812. Tomado de Wikipedia

La obra de Minard comprende otros mapas bastante interesantes y que demuestran su capacidad creativa. Un gran mérito de Minard es haber fabricado sus mapas en una época sin computadores y sin herramientas para capturar y procesar información geográfica. Actualmente existen varias opciones para hacer mapas y tal vez la más popular sea utilizar un software de Sistemas de Información Geográfico (SIG) como QGIS. Este camino es bastante eficiente y permite realizar productos de gran calidad. Sin embargo, hay algo que personalmente encuentro molesto cuando utilizo SIGs y es que no está articulado de manera natural con todo mi flujo de trabajo. Dado que todo el procesamiento y análisis de datos lo hago en R, hacer un mapa en un SIG implica exportar los datos desde R y volverlos a cargar a QGIS. Adicional a esto, y dependiendo del software y tipo de mapa, puede llegar a ser muy complicado hacer mapas que sean reproducibles y que se puedan reciclar para hacer mapas similares.

Por supuesto, existen también varias opciones para hacer los mapas directamente en R como ggplot o la función plot(). El año pasado empecé a utilizar la librería tmap() desarrollada por Martijn Tennekes y en mi opinión es bastante fácil de usar, ofrece mucha flexibilidad y las funciones que trae permiten configurar muchos detalles de los mapas. En este blog voy a hacer una presentación de la librería y voy a enseñar las principales funcionalidades que he utilizado durante el último año. Espero que sea de gran ayuda y que con algo de esfuerzo sea más fácil producir mapas que, como el de la Marcha de Napoleón de Minard, narren historias y presenten información compleja en una imagen.

Hay muchos otros recursos (mirar la última sección) que enseñan a usar tmap. Mis dos principales motivaciones para hacer otro tutorial son que, por un lado, no existen muchos documentos en español que personas en latinoamérica que no leen en inglés pueden consumir. Y por el otro lado, quería hacer un documento para mi mismo con las principales operaciones que realizo en tmap y tener un guía personalizada para consultar en el futuro.

Introducción y Rudimentos

La lógica de tmap es similar a la de ggplot por lo que algo de experiencia en ggplot hará que sea mucho más fácil entender. Sin embargo, no es necesario tener conocimientos previos. De manera muy general, primero se define el objeto espacial que se va dibujar mediante la función tm_shape() y después se indica que se quiere graficar con ese objeto espacial y se pasa otra función dependiendo de su clase. Por ejemplo, si es un objeo de líneas la función a utilizar es tm_lines() o si es un objeto de polígonos se usa la función tm_polygons(). A continuacón se pueden agregar nuevos objetos espaciales usando de nuevo tm_shape() y usando una función para dibujarlo. Al agregar un nuevo objeto espacial este se dibuja sobre lo que ya se ha graficado, un funcionamiento análogo al de los SIGs.

Otra aspecto a tener en cuenta es que las funciones se van agregando mediante el operador “+”. Todas las funciones en tmap comienzan de la forma tm_.

Además de tmap, vamos a utilizar otras dos librerías:

library(tmap)
library(sf)

Primero creemos los objetos AlquilerBicicletas y Londres los cuales contienen puntos donde se puden rentar bicicletas en Londres y la división geográfica de Londres respectivamente. 

AlquilerBicicletas <- spData::cycle_hire
Londres <- spData::lnd

Las estaciones y la ciudad se pueden dibujar de la siguiente manera:

tm_shape(shp = AlquilerBicicletas) + tm_dots()

tm_shape(Londres) + tm_polygons()

El primer mapa parece una nube de puntos y el segundo se limita a una división administrativa de la ciudad. Sobreponiendo las estaciones en los polígonos y configurando algunos argumentos se obtiene un mapa más interesante. 

tm_shape(Londres) + tm_polygons(col = "ONS_INNER") +
  tm_shape(shp = AlquilerBicicletas) + tm_dots(size = "nbikes", alpha = 0.2, col = "blue")

En este nuevo mapa hicimos tres operaciones adicionales:

Esta representación permite ver que las estaciones están en el centro de Londres y sugiere que hay unos puntos de concentración de préstamos (hotspots). Casi en el centro del mapa hay una zona en la que al parecer no hay estaciones ¿alguna idea de por qué?

Un objeto de líneas se grafica siguiendo la misma lógica

Sena <- spData::seine #Rio Sena (y otros) en Francia
tm_shape(Sena) + tm_lines()

Modo Interactivo

Para hacer mapas interactivos hay que cambiar el modo plot al modo view en la función tmap_mode().

tmap_mode("view")

tm_shape(Londres) + tm_polygons(col = "ONS_INNER", alpha = 0.5) +
  tm_shape(shp = AlquilerBicicletas) + tm_dots(size = "nbikes", alpha = 0.2, col = "blue")
tmap_mode("plot")

Si el mapa no se muestra en el navegador puede correr el código de manera local. Al acercarse (zoom in) al área donde hay pocas estaciones se puede obsevar que hay varios parques.

Haciendo un mapa


Ya vimos la lógica y algunos aspectos de tmap, y ya estamos en la capacidad de hacer algunos mapas. Pero todavía faltan varios elementos para que nuestros mapas se vean mejor. Nos interesa configurar la leyenda que se muestra, agregar una brújula y mostrar la escala. Todo esto lo vamos a ilustrar usando datos del Síndrome de Muertes Infantíles Repentinas en Carolina del Norte. 

nc = st_read(system.file("shape/nc.shp", package="sf")) #st_read hace parte de sf y se está usando para leer un archivo shape. 
## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.5/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID):    4267
## proj4string:    +proj=longlat +datum=NAD27 +no_defs

La leyenda siempre aparece por defecto cuando configuramos algún argumento como una variable de algún objeto en tm_shape(). La brújula y la escala se crean usando tm_compass() y tm_scale_bar(). 

tm_shape(nc) + tm_polygons(col="BIR79") + 
  tm_compass() + tm_scale_bar()

leyendas - función tm_layout()

Mediante las funcionesn tm_layout() o tm_legend() se pueden hacer cosas como establecer el tamaño usando el argumento scale, fijar un marco con legend.frame = T o poner un título con title. 

tm_shape(nc) + tm_polygons(col="BIR79") + 
    tm_layout(scale=0.8,
              legend.frame = T,
              legend.position = c(0.2,0.03))

La posición de la leyenda en el mapa es un vector con dos elementos para las coordenadas en el eje X y en el eje Y. En el mapa anterior las coordenadas eran legend.position = c(0.2,0.03) y para el siguiente la uicamos en la esquina superior izquierda haciendo legend.position = c(0.005,0.69). También es importante notar que el titulo de la leyenda no se establece directamente en tm_layout() sino que se toma del argumento title en tm_polygons(). De esta manera, hay una subleyenda para cada variable que se utiliza dentro del mapa para expresar color, tamaño o alguna otra característica. 

tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84") + 
    tm_legend(legend.position = c(0.005,0.69),
              scale=0.61,
              legend.frame = T)

A la leyenda se puede agregar información que no necesariamente se ha usado dentro del mapa. En el mapa que se muestra a continuación se incluyeron dos subleyendas que representan áreas y puntos que ni siquiera existen en el mapa. La función tm_layout() tiene más de 30 argumentos por lo que solamente hacer ajustes en la leyenda ya proporciona muchas opciones y flexibilidad. 

tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84") + 
    tm_layout(title = "Este es mi titulo",
              legend.position = c(0.01,0.05),
              scale = 0.9,
              legend.frame = T) + 
   tm_add_legend(type = c("fill"),
                labels = c("Algún Área","Otra Área"), 
                col = c("#66c2a5","#e5c494"),
                title = "Otras áreas") +
  tm_add_legend(type = "symbol", 
                labels = c("Un punto","Otro punto"), 
                col = c("#e78ac3","grey"))

¿Dónde está el sur?

Pese a que por lo general en un mapa el norte siempre va a estar arriba, el sur abajo, el oriente a la derecha y el occidente a la izquierda; se acostumbra poner la brújula para referenciar los puntos cardinales. Saber configurar la brújula no solo permite seguir el estándar que a veces hace que los mapas se vean más profesionales, sino que es una herramienta indispensable cuando aparece la excepción a la regla y el norte ya no está en la parte de arriba. Esto puede pasar cuando las geometrías que se están dibujando no se ajustan muy bien al rectángulo del lienzo y se rotan los objetos espaciales, un caso frecuente es el del mapa de Bogotá (ver acá, acá o acá).

tm_compass() también tiene el argumento position y se configura exactamente igual que en tm_layout(). El argumento type permite escoger el tipo de brújula y size corresponde al tamaño. 

tm_shape(nc) + tm_polygons(col="BIR79") + 
    tm_layout(legend.position = c(0.01,0.65),
              scale=0.65,
              legend.frame = T) + 
    tm_compass(size = 5, type = "8star")

El argumento show.labels marca cuántas de las direcciones cardinales van a aparecer en la brújula y cardinal.directions permite modificar las direcciones. 

tm_shape(nc) + tm_polygons(col="BIR79") + 
    tm_layout(legend.position = c(0.01,0.65),
              scale=0.65,
              legend.frame = T) + 
    tm_compass(size = 5, type = "rose", lwd = 3, show.labels = 0) 

tm_shape(nc) + tm_polygons(col="BIR79") + 
    tm_layout(legend.position = c(0.01,0.65),
              scale=0.65,
              legend.frame = T) + 
    tm_compass(size = 5, type = "arrow", lwd = 3, show.labels = 2, cardinal.directions = c("S", "W", "N", "E"))

La Escala

tm_scale_bar() es muy parecida a las dos últimas funciones y tiene argumentos como breaks con los que se puede experimentar. 

tm_shape(nc) + tm_polygons(col="BIR79", title = "Births 1979 - 84") + 
    tm_layout(legend.position = c(0.01,0.65),
              scale=0.65,
              legend.frame = T) + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 2, position = c(0.06,0))

tm_shape(nc) + tm_polygons(col="BIR79", title = "Births 1979 - 84") + 
    tm_layout(legend.position = c(0.01,0.65),
              scale=0.65,
              legend.frame = T) + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 1, position = c(0.06,0),
                 lwd = 5, color.dark = "darkgreen", color.light = "blue")

Si quiere probar diferentes opciones es bueno revisar la documentatción para la leyenda, la documentatción para la brújula y la documentatción para la escala.

Por el momento ya podemos hacer un mapa de gran calidad usando código como el siguiente. Esta vez no se graficó el mapa de manera inmediata sino que se guardó como un objeto de clase tmap que se puede llamar más adelante o al que se le puden agregar nuevas funciones. 

MyMap <- tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84") + 
    tm_layout(legend.position = c(0.005,0.69),
              scale=0.61,
              legend.frame = T,
              bg.color = "lightblue") + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 1, position = c(0.06,0),
                 lwd = 3)
MyMap

Colores

Personalmente los colores siempre son la parte más difícil de definir para un mapa. Cada persona tiene sus gustos y cuando el mapa tiene diferentes capas es complicado encontrar una combinación de colores que mantenga la consistencia. La forma más natural es configrar el argumento palette de varias de las funciones que dibujan geometrías. 

tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84", palette = "Blues") + 
    tm_layout(legend.position = c(0.005,0.69),
              scale=0.61,
              legend.frame = T,
              bg.color = "lightblue") + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 1, position = c(0.06,0),
                 lwd = 3)

En la función palette se puede ingresar cualquiera de las paletes de R o las proporcionadas por RColorBrewer. También se pueden invertir agregando un símbolo negativo al inicio de la paleta (palette = “-YlGnBu”) 

library(RColorBrewer)
#display.brewer.all()
tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84", palette = "YlGnBu") + 
    tm_layout(legend.position = c(0.005,0.69),
              scale=0.61,
              legend.frame = T,
              bg.color = "lightblue") + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 1, position = c(0.06,0),
                 lwd = 3)

Paleta personalizada

Miremos bien cómo es el comportamiento de la variable que estamos usando para colorear los polígonos.

summary(nc$BIR79)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     319    1336    2636    4224    4889   30757
library(ggplot2)
ggplot(nc,aes(x=BIR79)) + geom_histogram()

Considerando que el comportamiento de la variables parece seguir una distribución Normal sesgada, los siguientes intervalos podrían tener más sentido. 

Intervalos <- c(0, 2000, 4000, 6000, 31000)
Etiquetas <- c("0 - 2000", "2000 - 4000", "4000 - 6000", ">6000")
tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84", palette = "YlGnBu",
                           breaks = Intervalos, labels = Etiquetas) + 
    tm_layout(legend.position = c(0.005,0.69),
              scale=0.61,
              legend.frame = T,
              bg.color = "lightblue") + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 1, position = c(0.06,0),
                 lwd = 3)

También se puede crear un objeto que tengo los colores que se quieren usar para cada intervalo y pasarlo en el argumento palette.

MyPalette <- c("#f2f0f7", "#cbc9e2", "#9e9ac8", "#6a51a3")
tm_shape(nc) + tm_polygons(col="BIR79", title = "Nacimientos 1979 - 84", palette = MyPalette,
                           breaks = Intervalos, labels = Etiquetas) + 
    tm_layout(legend.position = c(0.005,0.69),
              scale=0.61,
              legend.frame = T,
              bg.color = "lightblue") + 
    tm_compass(size = 5, type = "rose", lwd = 3) + 
    tm_scale_bar(breaks = c(0,50,100,150,200),size = 1, position = c(0.06,0),
                 lwd = 3)

Mapas Auxiliares

Un mapa auxiliar es un mapa que se ingresa dentro de otra mapa para complementarlo o dar contexto de ubicación. Por ejemplo, al mapa que hicimos de las estaciones de préstamos de bicicletas se le puede agregar un mapa de Reino Unido para que cuando las personas vean el mapa puedan deducir que estamos mostrando a Londres.

Creemos el Mapa de Reino Unido (UK).

Mundo <- spData::world
UK <- Mundo[Mundo$name_long=="United Kingdom",]
tm_shape(UK) + tm_polygons()

Para crear el mapa auxiliar hay que crear los mapas como objetos tmap y usar la función print() de la librería grid. 

MapaLondres <- tm_shape(Londres) + tm_polygons(col = "ONS_INNER") +
  tm_shape(shp = AlquilerBicicletas) + tm_dots(size = "nbikes", alpha = 0.2, col = "blue")

MapaUK <- tm_shape(UK) + tm_polygons()
library(grid)
MapaLondres
print(MapaUK, vp = viewport(0.17, 0.8, width = 0.25, height = 0.25))

Y por último, miremos la forma de ubicar el rectángulo de Londres en el mapa auxiliar.

CajaLondres <- st_bbox(Londres) %>% st_as_sfc()
tm_shape(CajaLondres) + tm_polygons(col = "white", border.col = "blue", lwd = 4)

#plot(st_geometry(MontevideoBox))
MapaUK <- MapaUK + tm_shape(CajaLondres) + 
  tm_polygons(col = "white", border.col = "blue", lwd = 1.5)
MapaLondres
print(MapaUK, vp = viewport(0.17, 0.8, width = 0.25, height = 0.25))

Temas Pendientes y otros recursos

Como dije en el comienzo, este documento es también una guía para mi mismo. Lo que implica que de manera deliberada y sesgada he dejado algunos temas por fuera. De todas maneras, creo que algunos de esos temas podrían ocupar por si mismos un documento independiente. Para las personas que quieran seguir explorando tmap les recomiendo revisar el Capítulo Octavo del libro de Geografía computacional con R. También está el get started! del creador de la librería, su artículo científico y por supuesto el repositorio en github donde pueden rastrear problemas y mirar los desarrollos nuevos que se están implementando.