miércoles, 8 de febrero de 2017

Cómo hacer un mapa de Chile en R

A continuación describo cómo construí un mapa de Chile en R, con el objetivo de mostrar cuáles eran las capitales regionales más caras y baratas de Chile para un ejecutivo ABC1.

El costo de cada ciudad se estimó mediante una medición de los precios de una canasta de productos y servicios consumidos por hogares ABC1, es decir, se trata de una canasta que toma en cuenta la calidad de los productos y servicios consumidos, independientemente de su precio. A partir de estos datos se construye un índice el ICVE (índice de costo de vida de ejecutivos), y se construyen grupos de ciudades de costo alto, costo intermedio y costo bajo (mediante un algoritmo de k-means clustering).

En el mapa interesa colorear las regiones según su clúster de ingreso, y añadir una etiqueta con el nombre de la capital regional y el costo de vida estimado.

Descargué un mapa vectorial de Chile del sitio web de la Biblioteca del Ciongreso Nacional:
http://www.bcn.cl/siit/mapas_vectoriales/index_html

En particular, me interesaba el mapa: Division regional: polígonos de las regiones de Chile.

Para construir el mapa utilicé los siguientes paquetes de R:

x <- c("sp", "RColorBrewer", "ggplot2", "ggmap", "rgdal", "rgeos", "maptools", "dplyr", "tidyr", "tmap")
lapply(x, library, lib.loc="C:/Users/Guillermo Acuña/Documents/R library", character.only = TRUE)


Luego cargué los datos del mapa vectorial y los datos que quería mostrar en el mapa:

#Cargar los datos del mapa vectorial
chile <- readOGR(dsn = "C:/Users/Guillermo Acuña/Google Drive/CEEN - Area Compartida/2 Productos regulares/ICVE/2017/Mapa Chile", layer = "division_regional")

#Cargar los datos que se quieren mostrar en el mapa
icve <- read.csv(file='C:/Users/Guillermo Acuña/Google Drive/CEEN - Area Compartida/2 Productos regulares/ICVE/2017/Mapa Chile/default.csv', header = TRUE)


Los datos a mostrar son los siguientes:



Luego fusioné ambos datos:

#Unir los datos que se quieren mostrar en el mapa con los datos del mapa vectorial
chile@data <- left_join(chile@data, icve)


Si se examina la base de datos (View(chile)), se observan algunas filas repetidas y algunas con valores NA. Por otro lado, quería separar el mapa en dos, uno que mostrara las 7 regiones más al norte de Chile, y otro que mostrara el resto, por lo que generé dos sub-conjuntos de la base de datos original:

chile1 <- chile[complete.cases(chile$Costo),]
chile1 <- chile1[complete.cases(chile1$COD_COM),]
chile1 <-  chile1[chile1$COD_REGI <= 7 | chile1$COD_REGI == 13 | chile1$COD_REGI == 15,]

chile2 <- chile[complete.cases(chile$Costo),]
chile2 <- chile2[complete.cases(chile2$COD_COM),]
chile2 <-  chile2[chile2$COD_REGI > 7 & chile2$COD_REGI != 13 & chile2$COD_REGI != 15,]
chile2$Ciudad[6] = NA
chile2$Ciudad[7] = NA


Luego generé una paleta de colores, para pintar el mapa:

Colores <- brewer.pal(3, "RdYlGn") # Tres colores, uno para cada nivel.

Después generé los gráficos (sin imprimirlos aún):

g3 = tm_shape(chile1) +
  tm_fill(title = "", col = "Costo", palette = Colores, legend.show = FALSE) +  
  tm_layout(outer.margins = c(0, 0, 0, 0), inner.margins = c(.0, .25, .0, .4), legend.outside = FALSE, frame = FALSE,
            legend.text.size = .5, legend.position = c("right","bottom"), legend.bg.color = "white",
            legend.bg.alpha = 1) +
  tm_borders(alpha=.5) +
  tm_text("Ciudad", size=.5, shadow=TRUE, bg.color="white", bg.alpha=0, alpha = 1, just = "left")


g4 = tm_shape(chile2) +
  tm_fill(title = "", col = "Costo", palette = Colores, legend.show = TRUE) +  
  tm_layout(outer.margins = c(0, 0, 0, 0), inner.margins = c(.0, .25, .0, .4), legend.outside = FALSE, frame = FALSE,
            legend.text.size = .5, legend.position = c("left","bottom"), legend.bg.color = "white",
            legend.bg.alpha = 1) +
  tm_borders(alpha=.5) +
  tm_text("Ciudad", size=.5, shadow=TRUE, bg.color="white", bg.alpha=0, alpha = 1, just = "left")


Para poder imprimir ambos gráficos en una sola imagen, utilicé una función que encontré en Google:

multiplot <- function(..., plotlist=NULL, cols) {
  require(grid)
 
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
 
  numPlots = length(plots)
 
  # Make the panel
  plotCols = cols                          # Number of columns of plots
  plotRows = ceiling(numPlots/plotCols) # Number of rows needed, calculated from # of cols
 
  # Set up the page
  grid.newpage()
  pushViewport(viewport(layout = grid.layout(plotRows, plotCols)))
  vplayout <- function(x, y)
    viewport(layout.pos.row = x, layout.pos.col = y)
 
  # Make each plot, in the correct location
  for (i in 1:numPlots) {
    curRow = ceiling(i/plotCols)
    curCol = (i-1) %% plotCols + 1
    print(plots[[i]], vp = vplayout(curRow, curCol ))
  }
 
}


Por último, imprimí el gráfico utilizando la función anterior:

multiplot(g3, g4, cols = 2)

El resultado es el siguiente: