23 May 2022 ¿Cómo diseñar un cartograma con R?

Un poco de contexto sobre R

Tiempo atrás, en el blog de UNIGIS publicamos un artículo sobre los cartogramas. El artículo mostraba cómo generar este tipo de representación cartográfica, utilizando un software libre como QGIS. En esta ocasión, y con dos sencillos conjuntos de datos, te vamos a mostrar cómo diseñar tus propios cartogramas con R. De todos modos, el cartograma final, no será más que una excusa para mostrar algunos paquetes y funciones básicas de R. ¿Quieres saber cómo diseñar un cartograma con R? ¡Vamos allá!

R es, para quien no haya oído hablar de él, un lenguaje de programación orientado a la computación estadística y al diseño de gráficos. Este lenguaje cuenta con algo más de dos millones de usuarios alrededor del mundo, según datos del propio consorcio. R está disponible como Software Libre bajo los términos de la Licencia Pública General GNU de la Fundación del Software Libre en forma de código fuente.

Claro vencedor en su disputa particular con otros lenguajes de programación orientados a la estadística como SAS o SPSS, de un tiempo a esta parte, su dominio se ha visto ligeramente ensombrecido por otro lenguaje de programación, como Python:

Un punto de partida para el diseño del cartograma

En ocasiones, tendremos la necesidad de representar un atributo de manera un tanto diferente al modo tradicional. Mas allá del típico mapa de coropletas, sea cual sea el método de clasificación aplicado, podemos querer emfatizar una determinada información, alterando la forma original de las cosas -que para el ejemplo escogido, será el tamaño y la forma de los países-. Y para llevar a cabo esta tarea, nada mejor que un cartograma.

Así pues, a lo largo de las siguientes líneas, veremos de forma pautada, cómo representar el número de artículos por país de origen, que responden a la búsqueda «QGIS – Archeology», en las primeras quince páginas de Google Scholar.

Para diseñar nuestro cartograma, precisamos pues de dos conjuntos de datos:

• Una capa de polígonos que contenga los países del mundo
• Un fichero que contenga la información a representar temáticamente.

Por lo que respecta a la base cartográfica a utilizar, vamos a trabajar con un fichero en formato GeoJSON, que se puede descargar desde este repositorio de github . También se utilizará un fichero en formato .xlsx, de elaboración propia, y que guarda los resultados de la búsqueda con Google Scholar.

Preparar el entorno de R

R, al tratarse de un lenguaje que se ejecuta mediante línea de comandos, se puede trabajar directamente desde la consola o terminal del ordenador en el cual esté instalado. De todos modos, lo más recomendable, es hacerlo desde RStudio: un entorno de desarrollo integrado para R. Este es un entorno multiplataforma que a buen seguro, se convertirá en un elemento indispensable, si te dispones a trabajar con R por vez primera.

Este artículo no pretende centrar el foco en RStudio, sino que asumiremos que éste se encuentra ya instalado en nuestro ordenador de trabajo. En RStudio, lo más práctico no es trabajar directamente sobre la consola, sino hacerlo a través de un script. En este script, se irán introduciendo todas las órdenes que vamos a ejecutar (o que podemos ir ejecutando sobre la marcha), con sus respectivos comentarios cuando sean necesarios. Así, será posible además reaprovechar una y otra vez este mismo código, modificando únicamente los datos de entrada.

Paso #1: La instalación y activación de los paquetes

En R, el flujo de trabajo básico se centra en instalar el paquete o paquetes que se van a necesitar a lo largo del proyecto. A continuación, se activan y se ejecutan las funciones que estos paquetes proveen. Pero, ¿Sabías que en 2020 R contaba ya con más de 16.000 paquetes? Y si bien es cierto, tampoco es preciso que te asustes (aún). Como casi siempre, para llevar a cabo la mayor parte de las tareas que vas a realizar con R en tu día a día, te bastará con conocer y dominar apenas unas decenas de ellos. A continuación, una muestra de los que se van a utilizar en esta ocasión:

# instalar los paquetes necessarios
install.packages("sf")
install.packages("cartogram")
install.packages("readxl")
install.package("dplyr")
install.packages("ggplot2")
install.packages("tidyr")
install.packages("tmap")

# activar las librerías instaladas
library(cartogram)
library(sf)
library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr)
library(tmap)

El paquete sf, ofrece funciones para la lectura datos geográficos, gestionar proyecciones, entre otros. El paquete readxl permite leer una hoja de cálculo en formato .xls o .xlsx. Los paquetes dplyr y tidyr se utilizan para la manipulación de objetos de tipo tidy (tablas). Los paquetes ggplot2 y tmap, ofrecen funciones de diseño de gráficos y mapas. Y el paquete cartogram, como es obvio, sirve para la creación de los cartogramas.

Paso #2: La lectura de los datos

Una vez instalados y activados los paquetes necesarios para realizar las tareas, pasaremos a leer los datos y a guardarlos como nuevos objetos dentro de R:

# establecer el directorio de trabajo (opcional, aunque recomendable)
setwd("/home/lluis/R/")

# importar las geometrías y el fichero excel
mundo <- st_read("paises.geojson")
articulos <- read_excel("QGIS_ARQUEOLOGIA_ARTICULOS.xlsx")

El resultado de ejecutar la función st_read() sobre el GeoJSON que contiene los países, además de crear un nuevo objeto llamado mundo, arrojará el siguiente resultado en la consola de Rstudio:

Reading layer `mundo' from data source `/home/lluis/R/paises.geojson' using driver `GeoJSON'
Simple feature collection with 180 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -85.60904 xmax: 180 ymax: 83.64513
Geodetic CRS:  WGS 84

Esta respuesta, nos ofrece información muy básica sobre el nuevo objeto: número de columnas o variables, tipo de geometría y sistema de referencia de coordenadas.

A continuación, procedemos a realizar un comprobación un poco más completa acerca de la clase de los dos nuevos objetos, los nombres de las columnas, la verificación del sistema de referencia de coordenadas de la capa que almacena el contorno de los países, … Para ello, podemos valernos de las funciones class(), names(), str(), head(), tail(), …

> class(mundo)
[1] "sf"         "data.frame"
> class(articulos)
[1] "tbl_df"     "tbl"        "data.frame"
> names(mundo)
[1] "id"       "name"     "geometry"
> names(articulos)
[1] "PAIS"           "INVESTIGADORES"
> st_crs(mundo)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

Alternativamente, también puede utilizarse la función View() para ver la totalidad de la tabla. Llegados a este punto, puede que nos interese también ver las geometrías del objeto mundo. Para ello, vamos a crear nuestro primer plot. Uno muy básico y simple, pero suficiente por el momento. Así pues, ejecutaremos la siguiente función:

plot(mundo$geometry)

Paso #3: La reproyección de un objeto de tipo sf

Una de las particularidades del paquete cartogram, es que precisa trabajar con datos proyectados. El objeto mundo de que disponemos (sf o simple feature), se encuentra en coordenadas geográficas. Así pues, tocará reproyectarlo por ejemplo, al sistema de coordenadas proyectado Pseudo-Web Mercator. Nos valdremos en este caso de la función st_transform().

mundo3857 <- st_transform(mundo, 3857)
plot(mundo3857$geometry, main = "Mapa del mundo (EPSG:3857")

Paso #4: El filtrado de observaciones

La reproyección modifica la apariencia de las geometrías. Sin embargo, hay algo en este último gráfico, que «molesta» a efectos de representación de nuestra variable. Y es que a causa de la reproyección, el continente antártico, ha aumentado considerablemente su tamaño. Sabemos de antemano, que no hay ningún artículo cuyo origen del autor o autora, sea la Antártida. Así pues, lo más sencillo en este caso concreto, será eliminar esta geometría del objeto mundo. Para ello, nos valdremos de la función filter():

mundo3857 <- mundo3857 %>%
  filter(!name %in% 'Antarctica')

Con este filtro, seleccionaremos todas aquellas observaciones del objeto mundo, cuyo valor en el campo [name] no sea ‘Antarctica’. Así pasamos de 180 observaciones, a 179.

Para comprobar visualmente el efecto del filtrado, basta con ejecutar nuevamente la función plot(mundo3857$geometry), y observar el resultado:

Paso #5: La vinculación de los datos recopilados

Llegados a este punto, deberemos asignar a cada observación del objeto mundo (país), el número de artículos encontrados sobre el uso de QGIS en arqueología. No será ninguna sorpresa descubrir que este proceso se llevará a cabo, mediante un enlace de tablas. El atributo común en ambos objetos, para realizar el join, será el nombre de los países. Así pues, nos valdremos de la función left_join() del siguiente modo:

mundo_articulos <- left_join(mundo3857, articulos, by=c('name'='PAIS'))

Una vez realizado el enlace, el aspecto de la tabla resultante será el siguiente:

Paso #6: La gestión de los valores NA (missing values)

El paquete cartogram, como otras muchas funciones y paquetes, no es muy amigo de la presencia de missing values o valores desconocidos en los datos. De hecho, de no hacer nada con ellos, el objeto resultante del diseño del cartograma, sólo conservará las geometrías de los países para los cuales dispone de algún valor numérico. Lo sencillo, sería asignar de forma masiva el valor cero a los registros NA. Pero cartogram, tampoco es demasiado amigo de los valores cero a la hora de calcular la deformación de las áreas. Así pues, asignaremos el valor 0.1 a todos los registros NA. Utilizaremos en este caso, la función replace_na():

mundo_articulos_noNA <- mundo_articulos %>% replace_na(list(INVESTIGADORES = 0.1))
View(mundo_articulos_noNA)

Paso #7: Ahora sí. Vamos a por el cartograma…

Siguiente paso. Teniendo la tabla convenientemente preparada , ejecutaremos la función cartogram_cont(). Bastará con indicar cuál es el objeto base así como el número de iteraciones a realizar. En términos generales, a mayor número de iteraciones, mayor será la deformación aplicada a las formas originales. ¿Probamos con 3 iteraciones?

articulos_cartograma <- cartogram_cont(mundo_articulos_noNA, "INVESTIGADORES", itermax = 3)

¿Curiosidad por ver el resultado obtenido? Espera. Aguanta sólo unas pocas líneas de código más 😉

Paso #8: Desandando el camino de los missing values

Sabemos a ciencia cierta, que los valores 0.1 asignados a los valores NA originales, son cualquier cosa menos verdaderos. Así, una vez creado el cartograma, procederemos a convertirlos nuevamente en valores cero. Para ello, creamos un nuevo vector, en el cual convertiremos todos los valores 0.1 presentes en la columna [INVESTIGADORES], a 0. Posteriomente, añadiremos este vector como una nueva columna al objeto mapa_final. En este caso, utilizaremos las funciones replace() y mutate().

dtemp<-replace(articulos_cartograma$INVESTIGADORES, articulos_cartograma$INVESTIGADORES == 0.1, 0)
mapa_final <- articulos_cartograma %>% mutate(INVESTIGADORES = dtemp)

Paso #9: Diseñando el mapa definitivo

Pues ahora sí. Ya es momento de darle color y contexto al cartograma que acabamos de calcular, y que lleva por nombre mapa_final.

a) El paquete tmap.

En este primer ejemplo que se muestra a continuación, se ha utilizando el paquete tmap para el diseño del mapa temático. Como puede comprobarse, son muchas las funciones, argumentos y posibilidades que ofrece este paquete y, en buena medida, dependerá del tipo de mapa a diseñar, la tipología de atributo a mostrar, etc. En cualquier caso, este es solamente un ejemplo del resultado final. En el siguiente enlace, encontraréis un nutrido conjunto de ejemplos e información, sobre el paquete en cuestión.

tm_shape(mapa_final) +
  tm_style("gray") +
  tm_layout(main.title = "Artículos: Uso de QGIS en arqueología",
            main.title.size = 1,
            title = "Recuento basado en Google Scholar (2022)", title.size = 0.8)+
  tm_polygons(col = "INVESTIGADORES", style = "cont",
              breaks= c(1,5,10,20,40), palette = "YlOrBr",
              title = "Artículos según origen") +
  tm_scale_bar(position = c("left", "bottom"), width = 0.2) +
  tm_compass(position = c("0.001", "0.07"), size = 1.5)+
  tm_credits("UNIGIS Girona, 2022")

b) El paquete ggplot2

Más allá del paquete tmap, existen en R otras opciones para el diseño de mapas. Otro ejemplo clásico, es el paquete ggplot2. A continuación, se muestran las líneas de código necesarias para confeccionar un sencillo mapa basado en el cartograma creado con anterioridad. Como sucedía con tmap, en el siguiente enlace vais a encontrar toda la información sobre este enorme paquete orientado a la creación de gráficos, entre los cuales, se encuentran los mapas.

ggplot(data = mapa_final) +
  theme_set(theme_gray()) +
  theme(legend.position = "bottom") +
  xlab("Longitud") + ylab("Latitud") +
  scale_fill_distiller(name = "Artículos según origen", palette = "Spectral",
                       limits = c(0,45), breaks = c(1,5,10,20,40)) +
  geom_sf(aes(fill = INVESTIGADORES)) +
  ggtitle(label = "Artículos: QGIS aplicados en arqueología",
          subtitle = "Según Google Scholar (2022)") 

Paso #10: ¿Qué mas puede ofrecer el paquete cartogram?

Pues lo cierto es que, además de la función cartogram_cont() que se acaba de ver, el paquete dispone de dos funciones más –cartogram_dorling() y cartogram_ncont()– que ofrecen dos nuevos modos de representar, la realidad mediante cartogramas. En un caso se trata de sustituir todas las formas originales de los países por círculos mientras que, en el segundo caso, se representan los países, como polígonos o áreas no contiguas.

cartograma_dorling <- cartogram_dorling(mundo_articulos_noNA, "INVESTIGADORES", itermax = 3)
cartograma_ncont <- cartogram_ncont(mundo_articulos_noNA, weight = "INVESTIGADORES", k = 6)