Post 3 - Bicis

Introducción

Hecho con R Markdown

Posts anteriores

Posts anteriores

Post 1

Post 2

Post 3 - Acá abajo

Consideraciones previas

Esta página es estática, no tiene ninguna parte “activa” más allá de los controles del mapa. Para eso hay que usar Shiny. Otro día.

Iba a hacer el mismo planteo que en el post anterior pero ya que estamos aprovechemos para limpiar. Desde el vamos andar por ahí con una base de 2,8 millones de registros y 19 columnas es un poco incómodo. De la base original vamos a redistribuir la información así:

• Un frame (viajes2019) con los viajes y la información mínima necesaria: fecha, id_usuario, duracion, estación origen y destino (2,8 millones de registros)
• Un frame (estaciones2019) con las estaciones y su ubicación: número de estación, latitud, longitud, nombre de estación (422 registros)
• Un frame con la información de usuarios (acá no lo voy a implementar pero sale igual que el de estaciones si quieren probar)

Y dejarlos guardados para levantarlos según sea necesario. Esto facilita mucho el trabajo y la rapidez al probar código sin tener que esperar largos minutos.

Mapa básico con leaflet

Estaciones

Cargamos las librerías y los dos frames que vamos a necesitar: el de viajes y el de estaciones. Eliminando información duplicada pasamos de un archivo de 648 megas a uno de 108 megas con los viajes más otro de 200k con la información de las estaciones.

library(tidyverse)
library(ggplot2)
library(leaflet)
library(htmltools)
library(htmlwidgets)
library(webshot)
library(rgdal)
library(magick)

viajes2019 <- read.csv("~/bicis/viajes2019.csv", stringsAsFactors=FALSE)
estaciones2019 <- read.csv("~/bicis/estaciones2019.csv", stringsAsFactors=FALSE)

Vemos la base de estaciones inmediatamente. Se puede pasar por cada estación y sale el nombre.

leaflet(estaciones2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
  addMarkers(lng = ~lon,
             lat = ~ lat,
             popup = ~htmlEscape(nombre))

Ok esto es bastante modesto. Vamos con lo que quería mostrar al principio, el mapa con la información.

Trabajo por zonas

Antes de lo gráfico veamos los datos necesarios. Reconstruimos la información de las zonas:

estaciones2019 <- estaciones2019 %>% mutate(
       ZonaLat=sprintf("%02d",trunc((((lat+34.66263)*82.2096)*2))),
       ZonaLong=sprintf("%02d",trunc((((lon+58.51142)*64.3416)*2))),
       Zona=paste(ZonaLong, ZonaLat, sep = "-"),
 )

Desarrollando acá tenemos que voy a dividir el área ocupada por las estaciones en una cuadrícula de 20x20. Para esto tomo los datos de la estación que está más al sur (-34.66263) y la que está más al este (-58.51142). Después tomo la diferencia entre esas y las que están más al norte/este, divido por 20 y listo. En el código está el número original y después multiplicado por dos por que el plan original era hacerlo en una cuadrícula de 10x10 pero quedaba demasiado grande todo. Estas “zonas” después se juntan y quedan ordenadas por pares tipo “05-18” o sea desde arriba y a la izquierda (cambia la referencia del mapa acá, esto para hacer más cómodas las cuentas) cinco cuadritos a la derecha y dieciocho para abajo. Es como una guía T. Cada división de la cuadrícula es de aproximadamente 800x800 metros.

Cuando hice circular el primer mapa tuve muchos comentarios respecto a que no se entendía. Por ahora lo dejo así por una cuestión de practicidad con el código. Van a ver que es muy fácil trabajar con zonas cuadradas para calcular las salidas y llegadas. Además agrupar las estaciones por zonas en vez de hacerlo a nivel individual o por comunas(que no tienen límites regulares) ayudan a comprender mejor donde se originan y hacia donde van los viajes.

Ahora bien al separar las bases perdimos la referencia inmediata de la zona. Hay que volver a unir las bases por estación así recuperan la información de la zona. Tenemos que unir las columnas “EstOrigen” y “EstDestino” con la de “num_estacion”. Como en este caso las columnas de los dos dataframes son de diferente nombre hay que especificarlo en el merge.

viajes2019 <- merge(viajes2019, estaciones2019, by.x=c("EstOrigen"),by.y=c("num_estacion"))

Vamos a tener que cambiarle los nombres a las columnas antes del próximo merge y sacar las que no vamos a usar. En este caso sólo nos quedamos con Zona.

viajes2019 <- viajes2019 %>% 
  rename(
    ZonaOri = Zona
    )

viajes2019 <- viajes2019 %>% select(-c(ZonaLat, ZonaLong, lat, lon, nombre))

Lo mismo para recuperar la zona de las estaciones de destino.

viajes2019 <- merge(viajes2019, estaciones2019, by.x=c("EstDestino"),by.y=c("num_estacion"))

viajes2019 <- viajes2019 %>% 
  rename(
    ZonaDes = Zona
    )

viajes2019 <- viajes2019 %>% select(-c(ZonaLat, ZonaLong, lat, lon, nombre))

Vemos los viajes por zonas.

zonaviaj2019 <- viajes2019 %>% select(fecha, ZonaOri, ZonaDes)
zonaviaj2019 <- zonaviaj2019 %>% mutate(ZonaViaj = paste(ZonaOri, ZonaDes, sep = "a"))
zonaviaj2019 <- zonaviaj2019 %>% group_by(fecha, ZonaViaj) %>% summarize(viajes = n())

zonaviaj2019 <- zonaviaj2019 %>% mutate(
                                 ZonaOriX = as.numeric(substr(ZonaViaj, 1, 2)),
                                 ZonaOriY = as.numeric(substr(ZonaViaj, 4, 5)),
                                 ZonaDesX = as.numeric(substr(ZonaViaj, 7, 8)),
                                 ZonaDesY = as.numeric(substr(ZonaViaj, 10, 11)),
                                 Dist = round(sqrt(abs(ZonaOriX - ZonaDesX) + abs(ZonaOriY - ZonaDesY)),1),
)

Distancias recorridas (aprox)

Con esto armamos una lista de viajes realizados ordenados por recorridos entre zonas. Vemos que muchos viajes se realizan dentro de una misma zona (distancia = 0). Podríamos graficarlo ya que estamos:

zonaviajdist2019 <- zonaviaj2019 %>% group_by(Dist) %>% summarize(viajtot = sum(viajes))
zonaviajdist2019 <- zonaviajdist2019 %>% mutate(distancia=trunc(Dist))
zonaviajdist2019 <- zonaviajdist2019 %>% group_by(distancia) %>% summarize(viajestotal = sum(viajtot))


ggplot(zonaviajdist2019, aes(x = distancia, y = viajestotal))+
scale_x_continuous(labels = scales::comma)+
scale_y_continuous(labels = scales::comma)+  
geom_bar(stat = "identity", fill = "orange")+
labs(title ="Viajes según distancia", x = "Distancia (unidades)", y = "Viajes")

Según esto la mayoría de los viajes se dan entre uno y dos casilleros de distancia. Más atrás los viajes dentro de la misma zona y por último los viajes que han tomado más distancia. No hay ninguno que sea de cinco casilleros o más lo que sugiere que no se usan las bicis para recorrer grandes distancias. Si bien no se cuenta con la información de recorridos específicos recordemos que aproximadamente cada zona es de 800x800 metros. Un movimiento de cinco casilleros puede ser, por ejemplo, un viaje en línea recta de tres kilómetros y medio. Existen otras formas de calcular esto, en un caso se hizo calculando la ruta más corta en bici mediante la API del mapa de la ciudad o con la información de google. Pero sin saber el recorrido exacto creo que es más representativa la aproximación por ahora.

Agregando características al mapa

Viajes por cada zona

Volviendo al mapa este sería el mapa básico contando todos los viajes al inicio (o sea al momento de tomar la bici). La lógica básica es ubicar cada estación en una zona para después agruparlas y contar cada zona como una unidad. Después se cuenta y se escala.

Generamos el frame con la información que necesitamos: viajes acumulados por cada zona, en este caso contamos los viajes en el origen. En general la operación para contarlos a la llegada es exactamente la misma.

zonaorigtotal2019 <- zonaviaj2019 %>% mutate(zonaorig = paste(sprintf("%02d",ZonaOriX), sprintf("%02d",ZonaOriY), sep = "-")) %>% group_by(zonaorig) %>% summarize(viajes = sum(viajes))

zonaorigtotal2019 <- zonaorigtotal2019 %>% mutate(
                                ZonaOriX = as.numeric(substr(zonaorig, 1, 2)),
                                ZonaOriY = as.numeric(substr(zonaorig, 4, 5)),

)

Y ahora hacemos el mapa usando las zonas para reconstruir la información de latitud/longitud. Primero definimos la paleta (pal) para poder aplicar la escala con colores. Le pasamos la información del rango a cubrir (en este caso la cantidad de viajes). En cuanto al mapa para que salga bien en el knit hay que especificar un par de opciones adicionales como recordarle el proveedor de tiles (lo que se ve en el fondo del mapa) y su ubicación. Finalmente se agrega la leyenda y se crean los rectángulos con la información del dataframe.

pal <- colorNumeric(
  palette = c("blue", "orange","red"),
  domain = zonaorigtotal2019$viajes)

leaflet(zonaorigtotal2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
   setView(-58.43, -34.6,zoom = 12)%>%
  addLegend("bottomright", pal = pal, values = ~zonaorigtotal2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    opacity = 1)%>%
  addTiles() %>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorigtotal2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorigtotal2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorigtotal2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorigtotal2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(viajes)
  )

Viajes con zona más indicación de viajes

El resultado es un mapa con la información encima. Pero podemos hacerlo más atractivo (en la versión interactiva). Comenzamos con un pop up que indique el nombre de la zona y cantidad de viajes.

leaflet(zonaorigtotal2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
   setView(-58.43, -34.6,zoom = 12)%>%
  addLegend("bottomright", pal = pal, values = ~zonaorigtotal2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    opacity = 1)%>%
  addTiles() %>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorigtotal2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorigtotal2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorigtotal2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorigtotal2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(viajes),
    highlightOptions = highlightOptions(color = "white", weight = 2,
      bringToFront = TRUE),
    popup = ~htmlEscape(paste("Zona:", zonaorig, "Viajes:", viajes, sep = " "))
  )

Mostrar estaciones

Como referencia sería útil marcar las ubicaciones de las estaciones así que también las agregamos al mapa.

leaflet(zonaorigtotal2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
   setView(-58.43, -34.6,zoom = 12)%>%
  addLegend("bottomright", pal = pal, values = ~zonaorigtotal2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    opacity = 1)%>%
  addTiles()%>%
      addMarkers(lng = ~estaciones2019$lon,
             lat = ~estaciones2019$lat,
             popup = ~htmlEscape(estaciones2019$nombre))%>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorigtotal2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorigtotal2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorigtotal2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorigtotal2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(viajes),
    highlightOptions = highlightOptions(color = "white", weight = 2,
      bringToFront = TRUE),
    popup = ~htmlEscape(paste("Zona:", zonaorig, "Viajes:", viajes, sep = " "))
  )

Mostrar estaciones pero más lindo

Está bien en teoría pero es imposible de usar bien con los markers (los indicadores de posición) estándar. Usaremos unos más apropiados y vamos a agrupar para que no se acumulen al alejar el mapa.

biciIcon <- makeIcon(
  iconUrl = "https://image.flaticon.com/icons/svg/424/424370.svg",
  iconWidth = 50, iconHeight = 50
 )

leaflet(zonaorigtotal2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
   setView(-58.43, -34.6,zoom = 12)%>%
  addLegend("bottomright", pal = pal, values = ~zonaorigtotal2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    opacity = 1)%>%
  addTiles()%>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorigtotal2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorigtotal2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorigtotal2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorigtotal2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(viajes),
    highlightOptions = highlightOptions(color = "white", weight = 2),
    popup = ~htmlEscape(paste("Zona:", zonaorig, "Viajes:", viajes, sep = " "))
  )%>%
      addMarkers(
            lng = ~estaciones2019$lon,
            lat = ~estaciones2019$lat,
            popup = ~htmlEscape(estaciones2019$nombre),
            icon = biciIcon,
            clusterOptions = markerClusterOptions()
            )

Mucho mejor. Como vamos a querer mostrar muchas cosas en este mapa (ok no muchas, dos o tres pero síganme igual) tenemos dos vías para seguir adelante: generar una animación para mostrarlo en capas o hacerlo dentro del mismo mapa. Pros&Cons: el mapa es más flexible pero no funciona en todas las páginas (hola Medium). El gif animado es limitado (sólo va a mostrar los mismos datos en secuencia) pero funciona en todas partes.

Mapas con control de capas

Mapa con opción de ver estaciones

Vamos con el mapa primero. Aclaro que si bien el mapa se ve “dinámico” en cuanto a que se puede mover, hacer zoom, explorar y todo eso no deja de ser un control estático. No se actualiza con datos nuevos y no tiene otra lógica más allá de la propia del mapa. Para generar aplicaciones que funcionen en la web se necesita otro tipo de framework como el que provee Shiny por ejemplo.

Sobre la base del mapa anterior agregamos algunos controles adicionales. El primero nos va a permitir elegir si mostrar o no las estaciones.

biciIcon <- makeIcon(
  iconUrl = "https://image.flaticon.com/icons/svg/424/424370.svg",
  iconWidth = 50, iconHeight = 50
 )

leaflet(zonaorigtotal2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
   setView(-58.43, -34.6,zoom = 12)%>%
  addLegend("bottomright", pal = pal, values = ~zonaorigtotal2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    opacity = 1)%>%
  addTiles()%>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorigtotal2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorigtotal2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorigtotal2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorigtotal2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(viajes),
    highlightOptions = highlightOptions(color = "white", weight = 2),
    popup = ~htmlEscape(paste("Zona:", zonaorig, "Viajes:", viajes, sep = " "))
  )%>%
      addMarkers(
            lng = ~estaciones2019$lon,
            lat = ~estaciones2019$lat,
            popup = ~htmlEscape(estaciones2019$nombre),
            icon = biciIcon,
            clusterOptions = markerClusterOptions(),
            group="Estaciones"
            )%>%
  addLayersControl(overlayGroups = c("Estaciones"), 
                   options = layersControlOptions(collapsed = FALSE))

Preparando los datos necesarios

Antes de seguir necesitamos volver y repasar los datos con los que estamos trabajando ahora. El dataframe zonaorigtotal2019 tiene la información ya agregada, o sea sumada por cada zona de forma que no podemos saber en que dia se realizaron estos viajes. Vamos a retomar el zonaorig2019 donde esta información ya está presente. Lo único que necesitamos ahora es separar de la columna fecha la información perteneciente al día de la semana y a la hora para después reagrupar según lo que queremos mostrar.

Por las limitaciones de las capas de leaflet (bah no son limitaciones sólo que estamos usando un control de capa para algo que lo excede un poco) lo único que vamos a poder generar es una opción excluyente más los adicionales. Entonces acá vamos a organizar el dataframe según si el día es hábil o no.

Esta línea

ZonaOri=paste(sprintf(“%02d”,ZonaOriX), sprintf(“%02d”,ZonaOriY), sep = “-”)

Convierte los datos de X e Y de la zona en la información de zona completa con un 0 adelante si el número es menor a diez. Si se hace con el número sin formatear la zona podría quedar así: “1-1” lo cual introduce errores más adelante ya que la lectura se realiza según la posición en la cadena y por lo tanto necesita que el formato sea “01-01”.

zonaorigdiahora2019 <- zonaviaj2019 %>%
  mutate(
    hora = as.numeric(substr(fecha, 4, 5)),
    diasemana = as.numeric(substr(fecha, 6, 6)),
    ZonaOri=paste(sprintf("%02d",ZonaOriX), sprintf("%02d",ZonaOriY), sep = "-")
  )

zonaorigdiahora2019 <- zonaorigdiahora2019 %>%
  mutate(
    habil = ifelse(diasemana>5,0,1)
)

zonaorigdiahabil2019 <- zonaorigdiahora2019 %>%
  group_by(ZonaOri, habil) %>%
  summarize(viajes = sum(viajes))

zonaorighabil2019 <- zonaorigdiahabil2019 %>% filter(habil == 1)
zonaorignohabil2019 <- zonaorigdiahabil2019 %>% filter(habil == 0)


zonaorighabil2019 <- zonaorighabil2019 %>% mutate(
                                ZonaOriX = as.numeric(substr(ZonaOri, 1, 2)),
                                ZonaOriY = as.numeric(substr(ZonaOri, 4, 5)),

)


zonaorignohabil2019 <- zonaorignohabil2019 %>% mutate(
                                ZonaOriX = as.numeric(substr(ZonaOri, 1, 2)),
                                ZonaOriY = as.numeric(substr(ZonaOri, 4, 5)),

)

Mapa con control de días hábiles/no hábiles, estaciones y bicisendas

Nada nos separa de nuestro mapa ahora. Lo único que tenemos que hacer es duplicar los rectángulos tomando los datos de sus respectivos dataframes y asignándolos a los grupos correspondientes. Ya que estamos agregamos una capa más con la información de bicisendas (descargada de https://data.buenosaires.gob.ar/dataset/ciclovias). Para esto necesitamos la librería rgdal así importamos el shapefile muy fácilmente.

biciIcon <- makeIcon(
  iconUrl = "https://image.flaticon.com/icons/svg/424/424370.svg",
  iconWidth = 50, iconHeight = 50
 )

bicisendas <- readOGR("~/bicis","ciclovias")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jose.ayala\Documents\bicis", layer: "ciclovias"
## with 2932 features
## It has 27 fields
## Integer64 fields read as strings:  id codigo alt_izqini alt_izqfin alt_derini alt_derfin cod_sent

bicisendas <- spTransform(bicisendas, CRS("+proj=longlat +ellps=GRS80"))


dominioviaj <- list(0, max(zonaorighabil2019$viajes, na.rm = FALSE))

pal <- colorNumeric(
  palette = c("blue", "orange","red"),
  domain = dominioviaj)

leaflet(zonaorighabil2019) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addTiles(urlTemplate = 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') %>%
   setView(-58.43, -34.6,zoom = 12)%>%
  addLegend("bottomright", pal = pal, values = ~zonaorighabil2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    group = "Días Hábiles",
    opacity = 1)%>%
  addTiles()%>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorighabil2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorighabil2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorighabil2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorighabil2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(zonaorighabil2019$viajes),
    popup = ~htmlEscape(paste("Zona:", zonaorighabil2019$ZonaOri, "Viajes:", zonaorighabil2019$viajes, sep = " ")),
    group="Días hábiles"
  )%>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorignohabil2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorignohabil2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorignohabil2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorignohabil2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(zonaorignohabil2019$viajes),
    popup = ~htmlEscape(paste("Zona:", zonaorignohabil2019$ZonaOri, "Viajes:", zonaorignohabil2019$viajes, sep = " ")),
    group="Días no hábiles"
  )%>%
  addPolygons(data = bicisendas, weight = 2, color = "green", group = "Mostrar bicisendas")%>%
      addMarkers(
            lng = ~estaciones2019$lon,
            lat = ~estaciones2019$lat,
            popup = ~htmlEscape(estaciones2019$nombre),
            icon = biciIcon,
            clusterOptions = markerClusterOptions(),
            group="Mostrar Estaciones"
            )%>%
            addLayersControl(baseGroups = c("Días hábiles","Días no hábiles"),
              overlayGroups = c("Mostrar Estaciones", "Mostrar bicisendas"), 
              options = layersControlOptions(collapsed = FALSE))

Mapas animados

El caso para la animación

Como decía más arriba si bien esto funciona los controles estándar de leaflet son limitados y si uno por ejemplo quiere mostrar que pasa a lo largo del día debería hacer 48 dataframes con este método, 24 para los días hábiles y otros tantos para los no hábiles. He hecho cosas malas con R pero todo tiene un límite así que pasemos a ver como podemos hacer una animación con esto.

El funcionamiento básico es iterar en un loop donde se agrupa el dataframe inicial según la variable que vayamos a trabajar, dibujamos el mapa en memoria, se graba en un archivo temporal (recuerden borrarlo después) y por último se retoman esos archivos en una animación.

Iterando para generar los mapas

Armamos el loop inicial para generar las imágenes de los mapas en nuestro directorio. Agregamos algunas etiquetas para ver el progreso.

zonaorighorahabil2019 <- zonaorigdiahora2019 %>% filter(habil == 1) %>% group_by(ZonaOri, hora) %>% summarize(viajes = sum(viajes)) 

dominioviaj <- list(0, max(zonaorighorahabil2019$viajes, na.rm = FALSE))

pal <- colorNumeric(
  palette = c("blue", "orange","red"),
  domain = dominioviaj)

for (i in 0:23) {
  zonaorighoratemp2019 <- zonaorigdiahora2019 %>% filter(hora == i & habil == 1) %>% group_by(ZonaOri, ZonaOriX, ZonaOriY) %>% summarize(viajes = sum(viajes)) 
  m <- leaflet(zonaorighoratemp2019, options = leafletOptions(zoomControl = FALSE)) %>%
    
  setView(-58.47, -34.58, zoom = 12) %>%
    
    addLabelOnlyMarkers(
    lng = -58.4, lat = -34.52,
    label = paste("Lunes a Viernes"),
    labelOptions = labelOptions(noHide = T, textOnly = T ,  textsize = "20px", style = list("font-style" = "bold")))%>%
addLabelOnlyMarkers(
    lng = -58.32, lat = -34.55,
    label = paste("Hora: ",i, sep = " "),
    labelOptions = labelOptions(noHide = T, textOnly = T ,  textsize = "20px", style = list("font-style" = "bold")))%>%
  addLegend("bottomright", pal = pal, values = ~zonaorighorahabil2019$viajes,
    title = "Viajes en bici",
    labFormat = labelFormat(prefix = ""),
    group = "Días Hábiles",
    opacity = 1)%>%
    addTiles() %>%
  addRectangles(
    stroke = FALSE,
    smoothFactor = 0, 
    fillOpacity = 0.6,
    lng1= -58.51142 + (zonaorighoratemp2019$ZonaOriX *(0.15542/20)), 
    lat1= -34.66263 + (zonaorighoratemp2019$ZonaOriY *(0.12164/20)),
    lng2= -58.51142 + ((zonaorighoratemp2019$ZonaOriX +1)*(0.15542/20)), 
    lat2= -34.66263 + ((zonaorighoratemp2019$ZonaOriY +1)*(0.12164/20)),    
    fillColor = ~pal(zonaorighoratemp2019$viajes)
  )
  
    saveWidget(m, 'temp.html', selfcontained = FALSE)
    webshot('temp.html', file=sprintf('habil%02d.png', i),
            cliprect = 'viewport')
}

Arreglo de las imágenes generadas

Como los mapas a veces quedan descentrados apelamos a la misma librería magick para procesar y eliminar lo que sobra. Barremos los archivos y hacemos un crop.

for (i in 0:23) {
  
  imagenhabil <- image_read(sprintf("habil%02d.png", i)) %>%  image_crop("692x604+300+140")
  image_write(imagenhabil, path = (sprintf("habilcor%02d.png", i)), format = "png" )
 
}

Creación de la animación

Finalmente tomamos los archivos que creamos en el último paso y hacemos la animación.

png.files <- sprintf("habilcor%02d.png", 0:23)

  
GIF.convert <- function(x, output = "mapa.gif")
        {
        image_read(x) %>%
        image_animate(fps = 1) %>%
        image_write(output)
        }

        GIF.convert(png.files)

Consideraciones finales

Que nos dice un mapa

El problema con los cuadrados es que en el papel (o en el R Studio mejor dicho) parecía una idea espectacular pero no sé si se entiende bien. Y la gracia de todo esto justamente es comunicar datos de forma sencilla, ahí capaz haya algo para mejorar. No obstante sí vale el ejercicio de haber encontrado algunas cuestiones como los patrones de uso que en un gráfico se ven bastante claras.

¿Por qué todo esto?

Quería aprender un poco de R y terminé dando una vuelta larga alrededor de varios conceptos. Creo que algo aprendí y la idea de todo esto es poder compartirlo por si a alguien le sirve. No es estrictamente un tutorial pero si un poco de código explicado (o al menos eso intenté).

El autor

Me llamo José Luis, soy sociólogo y analista de sistemas. Me interesan los campos de investigación de tecnología, datos y transporte aunque si aparece otro tema y me engancho nunca se bien donde puedo terminar.

Si llegaron hasta acá gracias por leer!