Como hacer un mapa muy bonito de España en ggplot2

(Este post apareció originalmente en mi blog).

Hace unas semanas leí un artículo en el cual Timo Grossenbacher mostraba como consiguió hacer, en mi opinión, uno de los mapas más bonitos que he visto nunca. Timo empleó la que es, en mi opinión, una de las librerias más expresivas y bellas que hay para hacer gráficos, ggplot2. La versión original de ggplot2 es para R, pero existe una versión de python no exhaustiva gracias a la buena gente de Yhat.

Asi que por supuesto, tenía que replicarlo.

Antes que nada, aquí está el mapa.

Mapilla

El código empleado para hacer el mapa lo podeis descargar en github. He compartido varias versiones del mapa para que se pueda observar como los diferentes cambios en las escalas afectan a la visualización.

Código.

Para empezar, importamos las librerías necesarias:

 
setwd("/DIRECTORIO_DE_TRABAJO/")

if (!require(rgdal)) {
install.packages("rgdal", repos = "http://cran.us.r-project.org")
require(rgdal)
}

if (!require(rgeos)) {
install.packages("rgeos", repos = "http://cran.us.r-project.org")
require(rgeos)
}
if (!require(rgdal)) {
install.packages("rgdal", repos = "http://cran.us.r-project.org")
require(rgdal)
}
if (!require(raster)) {
install.packages("raster", repos = "http://cran.us.r-project.org")
require(raster)
}
if(!require(ggplot2)) {
install.packages("ggplot2", repos="http://cloud.r-project.org")
require(ggplot2)
}
if(!require(viridis)) {
install.packages("viridis", repos="http://cloud.r-project.org")
require(viridis)
}
if(!require(dplyr)) {
install.packages("dplyr", repos = "https://cloud.r-project.org/")
require(dplyr)
}
if(!require(gtable)) {
install.packages("gtable", repos = "https://cloud.r-project.org/")
require(gtable)
}
if(!require(grid)) {
install.packages("grid", repos = "https://cloud.r-project.org/")
require(grid)
}
if(!require(tidyr)) {
install.packages("tidyr", repos = "https://cloud.r-project.org/")
require(tidyr)
}
}
 

El siguiente paso es importar los datos. Tras mucho buscar, encontré un archivo shapefile con los municipios españoles en ArcGis, sin ninguna atribución que pudiera encontrar.

Para obtener los datos del censo español, hice uso de la *fantástica* herramienta de extracción de datos proporcionada por el Instituto Nacional de Estadística. La herramienta es una pesadilla en términos de usabilidad, así que si queréis simplemente hacer el mapa he compartido los datos en el repositorio.

 
data_spain data_spain$municipality_code data_spain$People data_spain$Average.age

#Cargamos el shapefile y lo convertimos en un dataframe
municipalities_spain map_data_fortified_spain % mutate(id = as.numeric(id))

#Ahora unimos los datos del censo con los datos geométricos usando municipality_code como clave
map_data_spain % left_join(data_spain, by = c("id" = "municipality_code")) %>% fill(Average.age)
rm(data_spain)
rm(map_data_fortified_spain)
rm(municipalities_spain)
 

Finalmente, el código para hacer el mapa en sí. Hay muchísima lógica en dicho código orientada a hacer el mapa más bonito, os recomiendo mirar el artículo original para ver la evolución de los parámetros del gráfico, en particular todo lo relativo a la escala de colores.

 
# Aquí hacemos que los saltos de la escala de edades sean más bonitos e informativos visualmente

# encontramos los extremos del rango de edad
minVal maxVal # calculamos los valores de las etiquetas de los rangos de edad
labels brks # Redondeamos los extremos del rango de edad
for(idx in 1:length(brks)){
labels }

labels # definimos una nueva variable con los datos de la escala de edad
map_data_spain$brks breaks = brks,
include.lowest = TRUE,
labels = labels)

brks_scale labels_scale

theme_map theme_minimal() +
theme(
text = element_text(family = "Ubuntu Regular", color = "#22211d"),
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
# panel.grid.minor = element_line(color = "#ebebe5", size = 0.2),
panel.grid.major = element_line(color = "#ebebe5", size = 0.2),
panel.grid.minor = element_blank(),
plot.background = element_rect(fill = "#f5f5f2", color = NA),
panel.background = element_rect(fill = "#f5f5f2", color = NA),
legend.background = element_rect(fill = "#f5f5f2", color = NA),
panel.border = element_blank(),
...
)
}

#Esta función simplemente extiende los extremos de la escala de edad para llegar al mínimo y el máximo
extendLegendWithExtremes p_grob legend legend_grobs # grab the first key of legend
legend_first_key legend_first_key$widths # modify its width and x properties to make it longer
legend_first_key$grobs[[1]]$width legend_first_key$grobs[[1]]$x

# último valor de la leyenda
legend_last_key legend_last_key$widths

legend_last_key$grobs[[1]]$width legend_last_key$grobs[[1]]$x

# cambiamos también la posición de la última etiqueta para que no se superponga a la anterior
legend_last_label legend_last_label$grobs[[1]]$x

# Insertamos el nuevo color de la leyenda en la leyenda combinada
legend_grobs$grobs[legend_grobs$layout$name == "key-3-1-1"][[1]] <-
legend_first_key$grobs[[1]]
legend_grobs$grobs[legend_grobs$layout$name == "key-3-6-1"][[1]] <-
legend_last_key$grobs[[1]]
legend_grobs$grobs[legend_grobs$layout$name == "label-5-6"][[1]] <-
legend_last_label$grobs[[1]]

# Ahora lo mismo para el valor mínimo de la leyenda
new_first_label new_first_label$label new_first_label$x new_first_label$hjust

legend_grobs new_first_label,
t = 6,
l = 2,
name = "label-5-0",
clip = "off")
legend$grobs[[1]]$grobs[1][[1]] p_grob$grobs[p_grob$layout$name == "guide-box"][[1]]

# se usa esta función para dibujar la escala
grid.newpage()
grid.draw(p_grob)
}

p geom_polygon(data = map_data_spain, aes(fill = brks,
x = long,
y = lat,
group = group)) +
# municipality outline
geom_path(data = map_data_spain, aes(x = long,
y = lat,
group = group),
color = "white", size = 0.1) +
coord_equal() +
theme_map() +
theme(
legend.position = c(0.7, 0.03),
legend.text.align = 0,
legend.background = element_rect(fill = alpha('white', 0.0)),
legend.text = element_text(size = 14, hjust = 0, color = "#4e4d47"),
legend.title = element_text(size = 20),
plot.title = element_text(size = 28, hjust = 0.8, color = "#4e4d47"),
plot.subtitle = element_text(size = 20, hjust = 0.8, face = "italic", color = "#4e4d47"),
plot.caption = element_text(size = 14, hjust = 0.95, color = "#4e4d47"),
plot.margin = unit(c(.5,.5,.2,.5), "cm"),
panel.border = element_blank()
) +
labs(x = NULL,
y = NULL,
title = "Spain's regional demographics",
subtitle = "Average age in Spanish municipalities, 2011",
caption = "Author: Manuel Garrido (@manugarri) Original Idea: Timo Grossenbacher (@grssnbchr), Geometries: ArcGis Data: INE, 2011;") +
scale_fill_manual(
values = rev(magma(8, alpha = 0.8)[2:7]),
breaks = rev(brks_scale),
name = "Average age",
drop = FALSE,
labels = labels_scale,
guide = guide_legend(
direction = "horizontal",
keyheight = unit(2, units = "mm"),
keywidth = unit(70/length(labels), units = "mm"),
title.position = 'top',
title.hjust = 0.5,
label.hjust = 1,
nrow = 1,
byrow = T,
reverse = T,
label.position = "bottom"
)
)
extendLegendWithExtremes(p)
 

Este código está diseñado muy cuidadosamente para exportar una imagen con un ancho de 2400 píxeles.

Dado que las islas Canarias estan muy alejadas de la península, una práctica común es desplazar las islas más cerca de España, esto lo he hecho en Gimp.

Notas.

- Yo sabía que España tenía un problema poblacional, pero ¡madre mía! El noroeste del pais parece un gran asilo. El mapa original de Suiza tenia una escala de edad en el rango 40-52 años, pero he tenido que expandirlo a 40-56 debido al envejecimiento poblacional español.

- Una vez más, me he dado cuenta de lo mal que está el movimiento Open Data en España:

  • La herramienta de extracción de datos del INE parece que se hizo en los años 90 (tanto en usabilidad como en velocidad).
  • La información del censo de España se actualiza, ¡cada 10 años!. Esto significa que este mapa está usando la información más actualizada que existe de la población española, y es de 2011. Estos datos deberían actualizarse más frecuentemente en un mundo en el que todo cambia más rápido.
  • Si vais al mapa del artículo original, vereis que su mapa tiene una capa de topografía muy bonita encima de los municipios. Timo usó una imagen ráster a escala 1:1.000.000 con información topográfica proporcionada por la oficina federal de topografía suiza.
    Yo me tiré un dia entero buscando algo similar para España, y según parece el Centro Nacional de Información Geográfica sólo proporciona mapas ráster en la escala de 25 y 50 metros lo que obliga a descargarte cientos de archivos de imagen y unirlos luego. De verdad que yo no tenía ganas de hacer eso para hacer un mapa a escala tan pequeña. Al final, me hice mi propia imágen topográfica raster tomando imágenes de internet que procesé en Gimp. No incluyo el relieve en el mapa por que, al contrario que con Suiza (donde no hay municipios en la región de los Alpes), en España mostrar el relieve ocultaría la información de los municipios por debajo.

- Aunque tengo que decir que la cantidad de esfuerzo y dedicación que ha llevado realizar este mapa es impresionante (una de las principales razones por las cuales yo quería replicarlo), creo que debería haber una manera mejor de hacer gráficas customizadas en ggplot2. El código actual solo funciona con una resolución específica en mente, y se necesitan un montón de pruebas hasta llegar a encontrar los tamaños de fuentes que hacen que todo encaje. Idealmente ggplot2 usaría un sistema de referencia variable (algo como los em en css) que cambiara en funcion del tamaño del mapa.

Eso es todo por hoy, muchas gracias por vuestra atención 🙂

Tutorial de OpenLayers usando IPython, Brython y brythonmagic

No, no habéis leído mal, hoy vamos a hablar de OpenLayers, una librería javascript para hacer 'mapping' en el cliente (navegador).

Como sabéis, hablamos principalmente de Python porque nos gusta y porque nos divierte y este tutorial de una librería javascript lo vamos a realizar usando una sintaxis pythónica y, mientras aprendemos el uso básico de OpenLayers, también veremos algo de sintaxis javascript.

Para poder seguir este tutorial necesitaréis tener una versión moderna de IPython notebook instalada e instalar brythonmagic del que ya hablamos por aquí hace poco. El tutorial lo podéis descargar desde aquí.

Finalmente, quería agradecer a Roger Veciana por ofrecerse a revisar el tutorial y aportar mejoras importantes al mismo. Gracias monstruo.

Saludos.

P.D.: El tutorial está en inglés, si alguien no es capaz de seguirlo que avise por aquí y si el tiempo lo permite intentaré hacer una versión en la lengua de Cervantes.

Usando el GPS sin perdernos (usando pyproj)

Como todos sabéis, la tierra no es una esfera perfecta ya que su forma es irregular y sus polos están achatados. Esto complica la vida en muchos ámbitos de la vida como, por ejemplo, conocer una posición sobre la superficie terrestre con el menor error posible.

Está bien, Houston, tenemos un problema. Vamos a ver una serie de conceptos básicos, vamos a plantear un problema y lo vamos a resolver con la ayuda de nuestro lenguaje de programación preferido.

[Para este artículo se ha usado python 2.7.2 y pyproj 1.9.2]

Conceptos básicos

La tierra no es esférica (ni plana, como algunos creen aun), ya lo sabemos (algunos). Tampoco tiene una forma regular por lo que no existe una única definición de la forma de la tierra. Dependiendo de la zona del mundo en que nos encontremos deberemos aproximar la forma de la tierra mediante un modelo que tenga los menores errores posibles en nuestra área de interés. El geoide sería la forma que tendría la Tierra si se midiera sobre el nivel del mar. Esto es, una superficie ondulada que varía unos 100 m arriba y debajo de una elipsoide bien ajustada. Las elevaciones y las líneas de contorno de la tierra están refereridas al geoide y no al elipsoide, mientras que la latitud y la longitud están refereridas al elipsoide. ¿Por qué es ondulada la forma del geoide? Debido a irregularidades en el campo gravitatorio.

Por otra parte, normalmente, no llevamos un globo terráqueo con nosotros ya que no nos cabe en un bolsillo, llevamos mapas, que son más cómodos de transportar para irnos de excursión al monte. ¿Dónde quiero llegar con esto? Para el que no se haya dado cuenta, acabamos de convertir un objeto tridimensional en un objeto bidimensional y esto provoca que haya distorsiones. Es imposible que al pasar de 3D a 2D seamos capaces de mantener las mismas distancias, direcciones, formas y áreas. Al menos uno de estos parámetros cambiará (se distorsionará) al plasmarlo en un mapa 2D. Y llegamos al punto en el que entran en juego las proyecciones geográficas. Proyecciones geográficas hay muchas, algunas conservan áreas, otras distancias,...

Jo, siento que lo he explicado fatal pero es que el tema no es sencillo de explicar en dos párrafos. Un esquema de lo que he querido decir en los anteriores párrafos sería el siguiente:

Problema

Bueno, veamos el problema a ver si se entiende mejor lo que quiero mostrar:

Imaginad que hemos subido al monte, está nevado y nos hemos perdido. Llevamos nuestro receptor GPS que nos da la información en coordenadas geográficas refereridas al sistema de coordenadas WGS84 (todos los GPS funcionan con WGS84). Queremos que nos rescaten pero sabemos que la persona que nos puede rescatar solo tiene mapas en ED50 en coordenadas cartesianas (proyección UTM). Como son datums diferentes puede que si le paso las coordenadas en WGS84 no me encuentre porque su mapa esté desplazado (ya que está en otro datum). ¿Cómo puedo transformar mis coordenadas a las suyas para que me encuentre y no me congele? Exacto, usando pyproj.

Solución

Imaginad que me encuentro en el punto con longitud y latitud (-3.368620º, 37.054883º) en WGS84:

[googlemaps https://maps.google.com/maps?q=37.054883,+-3.368620&hl=en&sll=-3.36862,37.054883&sspn=0.016023,0.01929&t=h&ie=UTF8&ll=37.054835,-3.368597&spn=0.005994,0.00912&z=16&output=embed&w=425&h=350]

Vamos a transformar esas coordenadas con pyproj. Pyproj es una pequeña biblioteca que permite acceder a la biblioteca de proyecciones cartográficas proj.4 escrita en C. Pyproj permite hacer transformaciones de coordenadas y calcular círculos máximos. Vamos a usarla para hacer una transformación de coordenadas (ya llegamos al código python):

import pyproj
## http://pyproj.googlecode.com/svn/trunk/docs/pyproj.Proj-class.html
## Creamos el primer sistema de coordenadas en WGS84
## epsg:4326, http://spatialreference.org/ref/epsg/4326/
p1 = pyproj.Proj(init = "epsg:4326")
## Creamos el segundo sistema de coordenadas en ED50 UTM huso horario 30
## epsg:23030, http://spatialreference.org/ref/epsg/23030/
p2 = pyproj.Proj(init = "epsg:23030")
## Mi posición en WGS84
lon = -3.368620
lat = 37.054883
## Transformamos del sistema de coordenadas p1 (WGS84 en grográficas)
## al sistema de coordenadas p1 (ED50, UTM, Z30)
x2, y2 = pyproj.transform(p1, p2, lon, lat)
print x2, y2

Si lo hemos hecho todo bien el resultado que obtendremos será:

467327.137963, 4101226.2785 (las unidades son metros).

Y estas coordenadas serán las correctas en el mapa ED50 en UTM de nuestro salvador para que nos pueda encontrar sin ningún tipo de error. Imaginad por un momento que os habéis equivocado y en lugar de pasarle las coordenadas en ED50 UTM se las pasáis en WGS84 UTM. El resultado sería:

## epsg:4326, http://spatialreference.org/ref/epsg/4326/
p1 = pyproj.Proj(init = "epsg:4326")
## epsg:32630, http://spatialreference.org/ref/epsg/32630/
p2 = pyproj.Proj(init = "epsg:32630")
## Mi posición en WGS84
lon = -3.368620
lat = 37.054883
## Transformamos del sistema de coordenadas p1 (WGS84 en grográficas)
## al sistema de coordenadas p1 (ED50, UTM, Z30)
x2, y2 = pyproj.transform(p1, p2, lon, lat)
print x2, y2

Cuyo resultado sería:

467225.167945, 4101024.2792 (las unidades son metros).

Que estaría a unos 226 m al suroeste en el mapa de mi salvador. Imaginad lo trágico del error en un día de ventisca en la montaña nevada...

Algo menos trágico podría ser un error como el siguiente :-):

Miscelánea

Bueno, no sé si se ha entendido muy bien la complejidad del tema, no sé si he cometido algún error gordo en alguna explicación de más arriba pero os he presentado pyproj, una pequeña biblioteca con la que poder jugar. También queda claro que detrás de unas pocas líneas de código en python hay mucho tiempo, recursos y personas detrás de muy diversos campos, geógrafos, físicos, matemáticos, programadores,... Me gustaría resaltar que gracias a todo el trabajo de esta gente podemos hacer de forma sencilla lo que hemos visto hoy en esta entrada. Todo ello no es gratis, aunque lo podamos usar libremente. 😉

Hasta la próxima.

P.D.: Como siempre, tenéis los comentarios para corregirme, criticarme (de forma constructiva), mejorar las explicaciones (entre nuestros seguidores de twitter hay algunos mapólogos, ¿no? 1, 2, 3)...

P.D.2: En python tambien existen bindings para gdal y ogr, una librería más compleja y completa que pyproj, con lo que podríamos resolver este problema de otra forma. Algún día caerá algo también usando gdal y ogr.