De R a la web: Caso Espacial.

En este ejercicio publicaremos en la web un mapa. Haremos una versión estática y otra interactiva.

Preparando el archivo de Mapa

Descargue de este link la carpeta con el mapa, guardelo en su computadora.
Descomprima la carpeta. Los archivos crean un mapa de formato shapefile. Un mapa shapefile es complicado leer desde la nube, pues es un conjunto de archivos.
Cambie el mapa a formato JSON. Para ello vaya al mapshaper; desde ahí, suba los archivos de la carpeta, simplifiquelo, y luego exportelo como topojson.
Suba el archivo topojson a Github.
Guarde el link de descarga del archivo subido.

Versión estática

Traiga el mapa desde GitHub, necesita tener instaladas las bibliotecas sp, geojsonio y rgdal:

library(sp)
library(geojsonio)
library(rgdal)

fromGit="https://github.com/escuela-alacip/dataViz/raw/master/data/seattle.json"

wazipMap <- rgdal::readOGR(fromGit,stringsAsFactors = FALSE)

Puede ver el mapa aquí:

plot(wazipMap)

Decide qué vas a graficar:

El mapa presenta un polígono por cada zipcode. La idea es colorear cada polígono. Normalmente la data que uno desea colorear está en otro archivo. Ese archivo está aquí:

fromGit2="https://github.com/escuela-alacip/dataViz/raw/master/data/contriWA_2016.RData"

# lo traemos asi:
load(file=url(fromGit2))

Tenemos:

str(contriWA_2016,width = 60, strict.width = 'cut')

Esta data informa cuánto han contribuido diferentes personas a candidaturas políticas. Está por persona, NO por ZIP, entonces haremos un cálculo por zip. Por ejemplo, el promedio contribuido por zip:

# similar and saving space:
columnToAggregate=contriWA_2016$amount

# new ROW:
newROW_UNIT=list(zipCode=contriWA_2016$contributor_zip)

WA_zip_contri=with(contriWA_2016, 
                   aggregate(columnToAggregate, 
                             by=newROW_UNIT, 
                             mean)) #mean per zip code 

names(WA_zip_contri)[2]='AVE_Amount'

# obteniendo:
head(WA_zip_contri)

Combine la información con el mapa:

Nuestro mapa tiene estas variables:

names(wazipMap)

El zip code de cada poligono está en ZCTA5CE10. Esa servirá para el merge:

layerContrib=merge(wazipMap,WA_zip_contri,by.x='ZCTA5CE10', by.y='zipCode',all.x=F)

Pensemos en este mapa resultante como una capa. Por lo tanto, creemos una capa base, que sólo sean los límites del mapa:

library(rmapshaper)
# This will make just a border of the state
baseMap <- ms_dissolve(wazipMap)

La base será:

plot(baseMap)

Sobre esta base graficaremos la variable. La variable la organizaremos en cinco intervalos.

Instalemos las librerías que nos permitan seleccionar el color y crear los intervalos:

library(RColorBrewer)
library(classInt)

Definamos la variable a graficar:

varToPLot=layerContrib$AVE_Amount

Crear la secuencia de colores usando una paleta:

numberOfClasses = 5

colorForScale='YlGnBu'

colors = brewer.pal(numberOfClasses, colorForScale)

Producir el gráfico. Notese que los intervalos se producen con un método de clasificación conocido como quantile classification:

library(tmap)

creditsText="EPSG:4326 - Proj=longlat - datum=WGS84"

base= tm_shape(baseMap,projection = "longlat") + tm_polygons(col = 'black')

layer1= base +  tm_shape(layerContrib) + 
                tm_polygons("AVE_Amount", 
                            style="quantile",n=5,
                            title="Contributions", # titulo para legenda
                            palette=colorForScale,
                            border.alpha = 0.2) 

fullMap= layer1 + tm_compass(position = c('left','TOP'),type = 'arrow') +
                  tm_scale_bar(position=c("RIGHT", "BOTTOM"),width = 0.2)+
                  tm_credits(creditsText, position=c("left", "bottom")) 

fullMap

Ajustando posiciones:

fullMap +  tm_layout(main.title = "Choropleth",
                     main.title.position = 'center',
                     legend.position = c('RIGHT','center'),
                                    #bottom,left,top,right
                     inner.margins=c(0.1,0,0.1,0.2)) # para ganar espacio

En vez de crear intervalos, podríamos usar los valores de la variable a graficar para definir ciertos grupos. Por ejemplo, podria identificar el decile más alto y el más bajo:

quantile(layerContrib$AVE_Amount, c(.1,.9))

Entonces, puedo graficar sólo los grupos que cumplan con una condición:

#filtros:
top10=quantile(layerContrib$AVE_Amount, c(.9))
bot10=quantile(layerContrib$AVE_Amount, c(.1))

#nuevos Mapas!
mapBot=layerContrib[layerContrib$AVE_Amount<=bot10,]
mapTop=layerContrib[layerContrib$AVE_Amount>=top10,]

Ahora, usemos esos nuevos mapas como capas:

###
legendText="Areas to watch"
shrinkLegend=0.4
title="Top and Botton Average Contribution to elections in WA (2009-2023)"
###

base= tm_shape(baseMap,projection = 'longlat') + tm_polygons(border.alpha = 0)

layer_1= base +  tm_shape(mapTop) + 
                tm_polygons(col = 'green',border.col = NULL) 

layer_1_2= layer_1 + tm_shape(mapBot) + 
                tm_polygons(col = 'red',border.col = NULL) 

fullMap= layer_1_2 + tm_compass(position = c('left','TOP'),type = 'arrow') +
                  tm_scale_bar(position=c("RIGHT", "BOTTOM"),width = 0.2)+
                  tm_credits(creditsText, position=c("left", "bottom"))

fullMap

Añadiendo leyenda:

fullMap_leg= fullMap + tm_add_legend(type="fill",
                                     labels=c('good','bad'),
                                     col=c('green','red'),
                                     border.col=NA,
                                     title='to watch')
fullMap_leg

Acomodando posiciones:

fullMap_leg + tm_layout(main.title = "Highlights",
                        main.title.position = 'center',
                        legend.position = c('RIGHT','center'),
                                    #bottom,left,top,right
                        inner.margins=c(0.1,0,0.1,0.2))

Versión interáctiva

Necesitaremos usar leflet.

library(leaflet)

# color segun valor
paletteFun=colorQuantile("YlGnBu", 
                         varToPLot,
                         n = 5)

# the base map
base= leaflet() %>% addTiles()

final = base %>% 
         addPolygons(data=layerContrib,
                     weight = 1, #anchura de brde
                     opacity =  1, # 0 transparencia total
                     fillOpacity = 0.5, # contraste de paleta
                     fillColor = ~paletteFun(AVE_Amount)) # coloreando

final

Añadiendo leyenda:

final %>% addLegend(data=layerContrib,
                    "bottomright",
                    pal = paletteFun, 
                    values = ~AVE_Amount,
                    title = "Contributions",
                    opacity = 1)

Los valores por defecto salen como porcentaje, podemos ajustar ello:

final %>% addLegend(data=layerContrib,"bottomright", 
                    pal = paletteFun, 
                    values = ~AVE_Amount,
                    title = "Contributions",
                    opacity = 1,
                    # changes:
                    labFormat = function(type="quantile", cuts, p) {
                        n = length(cuts) # cuantas
                        lower=round(cuts[-n],2) # intervalos
                        upper=round(cuts[-1],2)
                        cuts = paste0(lower, " - ", upper) # nuevos limites
                        }
     )

De igual manera, podemos usar los nuevos mapas para solo colorear grupos particulares:

base= leaflet() %>% addProviderTiles("CartoDB.Positron") 
layer1= base %>%
        addPolygons(data=mapBot,
                    color='blue',
                    fillOpacity = 1,
                    stroke = F, # sin borde del poligono
                    group = "Bottom") # nombre

layer_1_2= layer1%>%addPolygons(data=mapTop,
                                color="red",
                                fillOpacity = 1,
                                group = "Top")

layer_1_2

A esta altura quizas te hayas dado cuenta que es dificil regresar al mapa inicial si has usado el zoom. Podemos añadir un botón para ello:

# indicando centro
textFun="function(btn, map){map.setView([47.751076, -120.740135], 7)}"

final= layer_1_2 %>%
    # adding the button
    addEasyButton(
        easyButton(icon="fa-home", # a symbol
                   title="Zoom to Level 1",
                   onClick=JS(textFun)))

final

Ahora podemos tener una leyenda interactiva:

final %>% addLayersControl(
        overlayGroups = c("Top", "Bottom"),
        options = layersControlOptions(collapsed = FALSE))

Guarde sus mapas en uno o más archivos en RPubs.

AUSPICIO:

El desarrollo de estos contenidos ha sido posible gracias al grant del Berkeley Initiative for Transparency in the Social Sciences (BITSS) at the Center for Effective Global Action (CEGA) at the University of California, Berkeley

RECONOCIMIENTO

EL Dr. Magallanes agradece a la Pontificia Universidad Católica del Perú, por su apoyo en la participación en la Escuela ALACIP.

El autor reconoce el apoyo que el eScience Institute de la Universidad de Washington le ha brindado desde el 2015 para desarrollar su investigación en Ciencia de Datos.

Curso: Visualización de Datos con enfoque reproducible de Data Science

De R a la web: Caso Espacial.