Clase 7: Mapas con Leaflet

Primero, instalamos los paquetes necesarios para esta sesión:

install.packages("leaflet")
install.packages("sp")
install.packages("chilemapas")

Cargamos los paquetes necesarios

Limpiamos el environment

rm(list=ls())

Mapas con Leaflet

Leaflet es una de las librerías abiertas más populares y reconocidas de JavaScript. El paquete de R con el mismo nombre nos permite crear mapas interactivos. Probablemente ya han visto mapas de leaflet antes, porque muchísimas páginas web las utilizan, como por ejemplo http://openstreetmap.org.

leaflet funciona por capas, y al igual que en ggplot existen muchísimas capas disponibles para mejorar el aspecto de nuestros mapas. Hoy veremos:

• Map Tiles: Elegimos el tipo de mapa que queremos. (addTiles, addProviderTiles)
• Markers: Marcadores específicos en un punto del mapa (latitud y longitud) (addMarkers (lat= lng=), addAwesomeMarkers)
• Polígonos: una región específica del mapa, como por ejemplo una comuna.

Para hacer nuestro primer mapa simple, vayamos a Google Maps para encontrar las coordenadas de la unviersidad. Aquí hacemos mapas con diferentes Provider Tiles

mapa <- leaflet() %>% 
  addProviderTiles(providers$Stamen.Toner) %>% 
  addAwesomeMarkers(lat=-33.486454, lng=-70.518628)

mapa

leaflet() %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addMarkers(lat=-33.486454, lng=-70.518628)

leaflet() %>% 
  addTiles() %>%   # Por defecto
  addMarkers(lat=-33.486454, lng=-70.518628)

Mapa de la RM con datos del MINEDUC

Ahora, vamos a aprender a agregar layers como marcadores o polígonos.

Abrir la base de datos Rendimiento_2018.csv. Esta representa el rendimiento escolar de alumnos de octavo básico de la Región Metropolitana, para el año 2018.

rend_2018 <- fread("Rendimiento_2018.csv")

Ahora podemos empezar a explorar nuestros datos. Primero hagamos un gráfico simple para ver la distribución de las notas en la muestra completa:

ggplot(rend_2018[prom_gral>=50 & prom_gral<=70], aes(x=prom_gral)) +
  geom_bar()

Podemos hacer un análisis más profundo y separar los datos por género. Con esto, podemos ver diferencias entre el rendimiento de niños y niñas. (Hombres gen_alu=1 y mujeres gen_alu=2)

ggplot(rend_2018[prom_gral>=50 & prom_gral<=70], aes(x=prom_gral)) +
  geom_bar(fill = "dark green") +
  labs(x="Promedio General", y="Número de alumnos", title="Distribución de notas por género", caption="Fuente: MINEDUC") +
  facet_wrap(facets = "gen_alu", ncol=1)

Ahora vamos a los mapas:

Primero, vamos a abrir los datos georeferenciales de las comunas de Santiago, desde el paquete de CRAN chilemapas. Este paquete tiene información para generar mapas de las divisiones políticas y adminsitrativas de Chile. Incluye información a nivel comunal, provincial y regional (entre otros). Es muy útil ya que nos entrega la geometría de cada unidad geográfica que quedamos mapear. En este caso, vamos a hacer un mapa de la RM por comuna. Una vez que cargamos el paquete, podemos crear objetos con los datos que necesitamos. En este caso, vamos a usar la función filter del paquete dplyr para quedarnos solo con las comunas de la Región Metropolitana. (si van a ??chilemapas encontrarán información de todos los datos que contiene este paquete)

comunas <- mapa_comunas %>% 
  filter(codigo_region==13)

Juntamos las comunas de los alumnos con las comunas del mapa.

rend_2018_geo <- merge(comunas, rend_2018, by.x="codigo_comuna", by.y="cod_com_alu")

Para identificar diferencias entre las comunas de Santiago, vamos a mapear el promedio general de los alumnos, por comuna (esta base de datos, por simplicidad, solo tiene a los alumnos que aprobaron). Esto nos mostrará una representación del promedio general promedio en cada comuna, en un mapa del Gran Santiago. Primero, creamos la variable de promedio general:

rend_2018_geo <- as.data.table(rend_2018_geo)
notas_geo <- rend_2018_geo[, nota:=mean(prom_gral), by=codigo_comuna]

Luego, con distinct nos quedamos solo con una observación por comuna, y seleccionamos solo las variables relevantes para crear el mapa en un nuevo objeto. Luego lo transformamos a formato sf para poder trabajar con leaflet

notas_geo <- distinct(notas_geo) %>% 
  select(codigo_comuna, geometry, nota) 
notas_geo <- st_sf(notas_geo)

Ahora, creemos un mapa. Con el paquete RColorBrewer podemos elegir muchas paletas de colores preexistentes. Con display.brewer.all() podemos verlas y elegir otra.

pal <- colorNumeric(    # definimos la paleta de colores
  palette = "YlOrBr",
  domain = rend_2018_geo$nota
)


mapa <- leaflet() %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addPolygons(data = notas_geo,
              fillColor = ~pal(nota),
              color = "#b2aeae",  # usamos código Hex para los colores
              fillOpacity = 0.7, # rellenamos con color los polígonos
              smoothFactor = 0.2,
              weight = 1) %>%   # para el grosor de la línea
  addLegend(pal = pal, # paleta de colores 
            values = notas_geo$nota,
            position = "bottomright",
            title = "Promedio General") %>% 
            addScaleBar(position = "topright") # Agregar escala arriba a la derecha

mapa

Vacunas contra COVID-19 en chile

Cargar la base de datos de las primeras dosis de las vacunas y eliminar las observaciones con comuna desconocida.

vacuna_dosis1 <- fread('https://raw.githubusercontent.com/MinCiencia/Datos-COVID19/master/output/producto80/vacunacion_comuna_1eraDosis.csv')
vacuna_dosis1<-vacuna_dosis1[!Comuna%like%"Desconocido",]

Debemos cambiar la estructura del dataframe para realizar el mapa, por ello, vamos a crear dos columnas, una de ellas serán las fechas de vacunación y la otra será la cantidad de personas vacunas en la fecha y comuna respectiva. Utilizaremos la función melt.

Columnas necesarias para hacer melt

Fechas <- names(vacuna_dosis1)
Fechas <- Fechas[6:124]

vacuna_dosis1_1<-melt(data = vacuna_dosis1,id.vars = c("Region","Comuna","Codigo comuna","Poblacion"),measure.vars = Fechas ,variable.name = "Fecha")

Porcentaje de primeras dosis por comuna

vacuna_dosis1_P <- vacuna_dosis1_1[,.(Suma_Primeras_dosis=sum(value)), by=.(`Codigo comuna`,Comuna,Poblacion)]
vacuna_dosis1_P[,Porcentaje_Primeras_dosis:=(Suma_Primeras_dosis/Poblacion),by=.(`Codigo comuna`,Comuna)]

Mapa comunas

Hagamos un mapa para ver el porcentaje de personas vacunadas por comunas.

vacuna_dosis1_P[, `Codigo comuna` := ifelse(`Codigo comuna` < 10000 , paste0("0",vacuna_dosis1_P$`Codigo comuna`), vacuna_dosis1_P$`Codigo comuna`)] # Agregar un cero adelante en codigo de las comunas que tengan 4 digitos 

mapa<-merge(mapa_comunas,vacuna_dosis1_P,by.x='codigo_comuna', by.y='Codigo comuna') # Merge entre el mapa de las comunas que nos proporciona el paquete chilemapas y el dataframe de los porcentajes. 

mapa<-st_sf(mapa) # Convertie el dataframe de data.table a sf

bins<-seq(0,1,0.2) # Intervalos del porcenateje para los colores del mapa

paleta<-colorBin(palette = 'Reds',domain = vacuna_dosis1_P$Porcentaje_Primeras_dosis,bins = bins) # Paleta de colores para el mapa

labels <- sprintf(
  "<strong>%s</strong><br/>%g Porc. Primera Dosis",
  mapa$Comuna, mapa$Porcentaje_Primeras_dosis
) %>% lapply(htmltools::HTML) # Labels al pasar arriba del mapa

leaflet(mapa)%>% 
  addProviderTiles(provider = providers$OpenStreetMap.Mapnik)%>%
  addPolygons(color = ~paleta(Porcentaje_Primeras_dosis),
              weight = 1, 
              fillOpacity = 0.8,
              label = labels)%>%
  addLegend(pal = paleta, values = ~bins, opacity = 1,position = "bottomright",title ="Porc. 1º dosis") # Agregar una leyenda

## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +no_defs).
## Need '+proj=longlat +datum=WGS84'

## Warning in paleta(Porcentaje_Primeras_dosis): Some values were outside the color
## scale and will be treated as NA

## Warning in paleta(Porcentaje_Primeras_dosis): Some values were outside the color
## scale and will be treated as NA