R: Mapas Temáticos (Choropleth)

Frecuentemente surge la necesidad de referenciar data con la ubicación geográfica. Dos insumos indispensables, para este propósito, son los datos de interés y las coordenadas geográficas (longitud, latitud) Con el paquete raster se pueden obtener rapidamente coordenadas geográficas, si el nivel de desagregación que se requiere no es muy detallado; es decir, si las regiones que se van a graficar involucran países o estados. La funcion raster::getData puede proveer data geográfica de cualquier parte, pero sólo en los niveles mencionados.
Por ejemplo, si el requerimiento fuera un mapa de suramérica, un modo de proceder sería el siguiente: obtenemos los nombres de los países que componen el subcontinente, la funcion getData requiere codificación ISO3 lo cual puede obtenerse mediante la función raster::ccodes() de la siguiente forma

library(sf)
library(gapminder)
library(ggplot2)
library(magrittr)
library(raster)

raster::ccodes()[ccodes()$continent == 'South America', 
                 c('NAME', 'ISO3')] %>% knitr::kable()

Para simplificar, uno puede excluir las islas.

c1 <- c('Uruguay', 'Brazil', 'Bolivia',
        'Colombia', 'Suriname', 'Ecuador',
        'Paraguay', 'Argentina', 'Chile',
        'Peru', 'Venezuela', "Guyana", "French Guiana")
        
cods <- ccodes()[ccodes()$NAME %in% c1, c('NAME', 'ISO3')]

cods %>% list2env(., envir = globalenv())

cods <- ISO3 %>%  as.list() %>% setNames(., NAME)

El objeto cods será una lista con nombres y códigos representando las regiones o países de interés, en este caso suramérica, luego para descargar las coordenadas se introduce este vector como argumento en getData

 coords_shp <- lapply(cods, 
function(x) getData(name = 'GADM', country = x, level = 0))

Una vez se tengan las coodenadas descargadas, pueden grabarse al disco duro con la función raster::shapefile:

Map(function(x, y) shapefile(x,
                             filename = paste(d1, 
                                              y, 
                                              sep = '/')), coords_shp, NAME)

rm(NAME, ISO3, cods)

Graficar el mapa con ggplot2 puede facilitarse leyendo los archivos (shapefiles) del disco con la ayuda de sf:

lst <- lapply(fnomb, function(x) read_sf(x)) # fnomb es la dirección de los archivos en el disco

lst <- lst %>% do.call(rbind, .) # lst se combinan en un único objeto
lst %<>% st_transform(., 29101) # proyección de las coordendas

Al asignar la importación de los archivos a lst el objeto se convierte en una lista, luego con do.call y rbind se convierten a un data.frame. Un modo de graficar el mapa, puede ser el siguiente:

lst %>% 
  ggplot() +  
  geom_sf(aes(geometry = geometry)) +
  geom_sf_text(aes(label = GID_0, geometry = geometry),
               fun.geometry =  st_centroid, size = 3) +
  theme_void()

Ahora solo falta encontrar algunos datos para referenciarlos en el mapa, gapminder ofrece algunos datos para países como 'Esperanza de Vida', 'PIB per capita', etc. A manera de ejemplo se presenta el Producto Interno Bruto per capita para el año 1997:

gp <- gapminder_unfiltered[gapminder_unfiltered$country %in% c1 & 
                       gapminder_unfiltered$year == 1997, ]
                       
gp[['gdpM']] <- gp$gdpPercap %>% 
  divide_by(., 10^3)%>% # convertir a miles
  base::cut(., breaks = nclass.Sturges(gp$gdpPercap)) # para convertir en categórica la variable
  
gp[['gdpMM']] <- paste(as.numeric( sub("\\((.+),.*", "\\1", gp$gdpM) ),
      as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", gp$gdpM) ),
      sep = "-") # remover los corchetes y parentesis

gp <- gp[, c('country', 'gdpMM')]

lst <- merge(lst, gp, by.x = 'NAME_0', by.y = 'country')

lst %>% 
  ggplot() +  
  geom_sf(aes(geometry = geometry, fill = gdpMM)) +
  geom_sf_text(aes(label = GID_0, geometry = geometry),
               fun.geometry =  st_centroid, size = 2,
               colour = "#CDAD00") +
  scale_fill_brewer(palette = 'Purples') +
  labs(caption = 'Gapminder: PIB per capita 1997\n($ miles)') +
  guides(fill = guide_legend(title = 'PIB\npercapita')) +
  theme_void()

Resulta un mapa igual al que esta al inicio de este post usando otros colores. La desventaja que encuentro en este procedimiento, es que el objeto que se obtiene al combinar todas las coordenadas es bastante grande, dependiendo del tipo de manipulación de la data que se tenga en mente puede ser improcedente. En todo caso, habrá otras opciones disponibles.