CIENCIA DE DATOS PARA LAS CIUDADES

Retyk, 2019

TP 1: Calculando y mapeando agregados por área + Adquiriendo open data urbana

Existe un debate en torno a la cantidad de m2 que deberían existir por habitante en la ciudad, y de cuan lejos estamos del celebre índice de la OMS de 10m2 de EEVV/hab.

Mas allá de la relación sup.EEVV/hab ¿Podemos ponderar cuan buenos son esos EEVV? Valorar la calidad de un Espacio Verde público es complejo porque se deben tomar en cuenta muchos componentes: balance de superficies, equipamiento, iluminación, cercanía a transporte público, cantidad de visitantes, asoleamiento, etc.

Para hacer una primera aproximación, vamos a tomar la cantidad de arboles en EEVV como indicador de calidad del espacio, partiendo de la simplificacion de que un EEVV con muchos árboles es mejor que uno con pocos.

Vamos a trabajar con los datasets de BA Datos Abiertos, usando el de árbolado en espacios verdes. Existe un dataset complementario de arbolado de alineacion (calles), pero por ahora usaremos solo el primero.

barrios <- st_read('http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios.geojson')
comunas <- st_read('https://bitsandbricks.github.io/data/CABA_comunas.geojson')
Arb.EEVV <- read.csv('C:/MEU (DiTella)/2019 - A2/2 Trimestre/MU115 - Ciencia de datos 2/Datos2/BA/Arbolado EEVV/arbolado-en-espacios-verdes.csv', encoding = "UTF-8")

Vamos a graficar la info de Árboles en EEVV

ggplot() +
  geom_sf(data = barrios) +
  geom_point(data = Arb.EEVV, 
               aes(x = long, y = lat),
               alpha = .3, 
               color = "Dark Olive Green")+
  
     labs(title = "Arbolado de Buenos Aires",
         subtitle = "Ciudad Autónoma de Buenos Aires, 2018")

Si queremos sectorizar los datos para analizarlos por cada barrio podemos realizar un join espacial.

Pero primero tenemos que filtrar el CSV para sacar las entradas sin coordenadas, en caso de que las hubiera, y luego transformamos el dataframe en uno espacial.

Arb.EEVV <- Arb.EEVV %>% 
    filter(!is.na(lat), !is.na(long)) %>% 
    st_as_sf(coords = c("long", "lat"), crs = 4326)

Realizamos el join espacial

EEVV.barrios <- st_join(Arb.EEVV, barrios)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Ahora que agregamos la información del barrio y comuna a cada registro de Arbol en espacios verdes, podemos analizar la cantidad de arboles en eevv por barrio o comuna.

Creamos un contador de Arboles en EEVV por Comuna

conteo.arbcom <- EEVV.barrios %>% 
    group_by(comuna, tipo_folla, origen) %>% 
    summarise(cantidad = n())

veamos las variables resultantes del nuevo dataset conteo

names(conteo.arbcom)
## [1] "comuna"     "tipo_folla" "origen"     "cantidad"   "geometry"

Lo vamos a usar para armar un grafico que muestre la cantidad por comuna, diferenciando especies exóticas de nativas.

La incorporación de especies nativas de la ecoregión a los espacios públicos es un debate que está tomando mucha relevancia en las ciudades, ya que tienen un menor consumo hídrico, se adaptan mejor a los cambios de temperatura y promueven la biodiversidad.

ggplot(conteo.arbcom) +
  
  geom_bar(aes(x = as.factor(comuna), weight = cantidad, fill = origen))+
      coord_flip()+
  
      labs(title = "Distribución de arboles en EEVV segun comuna",
       subtitle = "Composición por Origen",
       x = "Comuna",
       y = "Cantidad de Árboles",
       caption = "Fuente: Elaboración propia en base a datos de BA Data")+
  theme_minimal()

A grandes rasgos podemos observar que las comunas 1, 8 y 14 están muy por encima de las demás en cantidad de árboles por EEVV. Tambien podemos ver que en en la mayoría de las comunas hay una fuerte prevalencia de especies exóticas.

Para graficar esta info en el mapa primero vamos a realizar un join para volcar la cantidad a cada comuna

conteo.arbcom2 <- EEVV.barrios %>% 
    group_by(comuna) %>% 
    summarise(cantidad = n())
conteo.arbcom2 <- rename(conteo.arbcom2, "comunas" = "comuna")

Realizamos el join

comunas.v2 <- st_join(comunas, conteo.arbcom2)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Y finalmente graficamos la cantidad de arboles en EEVV por comunas

ggplot() +
    geom_sf(data = comunas.v2, aes(fill = cantidad)) +
    scale_fill_viridis_c()+
    
   geom_sf_label(data = comunas, 
                 aes(label = comunas)) +
  
   theme_void()+
 
     labs(title = "Cantidad de arboles en EEVV por Comunas",
         subtitle = "Ciudad Autónoma de Buenos Aires, 2018",
         fill = "Cantidad",
         caption = "Fuente: BA Datos")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
## not give correct results for longitude/latitude data

Hay que tener en cuenta que las comunas que mejor performan son las que tienen reservas ecológicas o grandes parques (Reserva Costanera sur, Reserva Lago Lugano, Parque Tres de Febrero). Las comunas 3 y 5, que son las que peor rankean, tienen muy pocos m2 de EEVV, y por lo tanto dificilmente puedan tener muchos árboles. Pero esto no nos habla del rendimiento de árboles por EEVV.

Para conseguir la relacion arbol/m2 verde, siguiendo el criterio que adoptamos al empezar el ejercicio, vamos a usar el dataset de espacios verdes.

EEVV <- read.csv('C:/MEU (DiTella)/2019 - A2/2 Trimestre/MU115 - Ciencia de datos 2/Datos2/BA/EEVV/espacio-verde-publico.csv', encoding = "UTF-8")

Obtenemos un dataframe resumen con los m2 verdes totales por comuna. Vamos a expresar los totales en hectáreas, para evitar trabajar con números tan grandes.

EEVV <- EEVV %>% 
    select("comunas", "area") %>%
    group_by(comunas) %>% 
    summarise(total_Ha_verde = sum(area) / 10000 )

Ordenamos por Superficie verde de mayor a menor

EEVV %>%
  arrange(desc(total_Ha_verde))
## # A tibble: 16 x 2
##    comunas total_Ha_verde
##      <dbl>          <dbl>
##  1       8        421.   
##  2       1        326.   
##  3      14        222.   
##  4      13        174.   
##  5      10        124.   
##  6      12        101.   
##  7       9         71.7  
##  8       4         61.3  
##  9       7         36.2  
## 10       2         33.1  
## 11       6         21.4  
## 12      15         20.8  
## 13      11         11.8  
## 14       3          3.48 
## 15       5          2.09 
## 16       0          0.718

Calculamos la superficie total

EEVV %>%
  summarise(Total = sum(total_Ha_verde))
## # A tibble: 1 x 1
##   Total
##   <dbl>
## 1 1632.

Graficamos la superficie verde por comuna

ggplot(EEVV) +
  
  geom_bar(aes(x = as.factor(comunas), weight = total_Ha_verde))+
      coord_flip()+
  
      labs(title = "Distribución de EEVV Público por comuna",
       subtitle = "Expresado en Hectáreas",
       x = "Comuna",
       y = "Superficie",
       caption = "Fuente: Elaboración propia en base a datos de BA Data")+
  theme_minimal()

Los resultados que arroja el dataset difieren de aquellos publicados por el Anuario Estadistico 2017 (Dirección General de Estadísticas y Censos del GCBA). Es posible que el anuario contemple todos los EEVV de la ciudad, y no solo los EEVV públicos (usados para este ejercicio).

Anuario 2017

Por Ultimo, vamos a unir la info de EV a la de arbolado

comunas.v3 <- left_join(comunas.v2, EEVV)
## Joining, by = "comunas"

Generamos el indice de Arboles Totales / Sup. Verde, para cada comuna

comunas.v4 <- mutate(comunas.v3, indice = cantidad / total_Ha_verde)

Graficamos

ggplot() +
    geom_sf(data = comunas.v4, aes(fill = indice)) +
    scale_fill_viridis_c()+
    
   geom_sf_label(data = comunas, 
                 aes(label = comunas)) +
  
   theme_void()+
 
     labs(title = "Indice Arbolado / Sup. Verde por Comunas",
         subtitle = "Ciudad Autónoma de Buenos Aires, 2018",
         fill = "Cantidad",
         caption = "Fuente: BA Datos")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
## not give correct results for longitude/latitude data

comunas.v5 <- comunas.v4 %>% st_set_geometry(NULL) %>%
  arrange(desc(indice))

comunas.v5$barrios <- as.factor(comunas.v5$barrios)
comunas.v5$barrios <- factor(comunas.v5$barrios, levels = comunas.v5$barrios[order(comunas.v5$indice)])
ggplot(comunas.v5) +
  geom_bar(aes(x=barrios, weight = indice, fill = indice)) +
  labs(title = "Indice Arbolado / Sup. Verde por Comunas",
       x = "Barrios",
       y = "Indice") +
  theme(axis.text.x = element_text(angle = 90, size = 6)) + 
  scale_fill_viridis_c()

Este enfoque evitaría penalizar a las comunas con poco espacio verde, y se centraría en la proporción arbolada de los EEVV, independientemente de la superficie de los mismos.

TP 2: Explorando y mapeando información georeferenciada de OpenStreetMap

El PROCREAR está terminando su proyecto más grande y ambicioso a nivel nacional, el conjunto de viviendas conocido como Estación Buenos Aires, en el límite entre Parque Patricios y Barracas. El complejo está conformado por 56 edificios que albergan 2369 viviendas de diferentes tipologías, que oscilan entre 30 y 100 m2, y se calcula que recibirán cerca de 10.000 nuevos vecinos.

Teniendo en cuenta las dimensiones del proyecto, y por la relevancia que tiene para la zona sur de la ciudad, seria oportuno relizar un estudio de la zona para

bbox <- getbb("Parque Patricios, Ciudad Autonoma de Buenos Aires")
bbox
##         min       max
## x -58.41340 -58.38989
## y -34.64862 -34.62723
bbox_poly <- getbb("Parque Patricios, Ciudad Autonoma de Buenos Aires", format_out = "sf_polygon")
leaflet(bbox_poly) %>%
    addTiles() %>% 
    addPolygons()
patricios <- opq(bbox) %>% 
    add_osm_feature(key = "highway")
patricios <- patricios %>% 
    osmdata_sf()
patricios
## Object of class 'osmdata' with:
##                  $bbox : -34.6486151,-58.4133954,-34.6272312,-58.3898933
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 1654 points
##             $osm_lines : 'sf' Simple Features Collection with 519 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 7 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL
calles <- patricios$osm_lines

head(calles)
ggplot() +
    geom_sf(data = calles)

Aplicamos el polygono para limpiar lo que esta por fuera

calles <- st_intersection(calles, bbox_poly)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

Y graficamos nuevamente

ggplot() +
    geom_sf(data = calles)

Ahora que tenemos las calles, vamos a usar la informacion que tienen

Vamos a usar maxspeed para ver la velocidad de las vias de circulación

Primero tenemos que pasar la info a dato numeral

calles <- calles %>% 
  mutate(maxspeed = as.numeric(maxspeed),
         lanes = ifelse(is.na(lanes), 1, as.numeric(lanes)))

Ahora graficamos

ggplot(calles) +
    geom_sf(aes(color = maxspeed), alpha = 0.5) +
    scale_color_viridis_c() +

    labs(title = "Parque Patricios",
         subtitle = "Vías de circulación",
         caption = "fuente: OpenStreetMap",
         color = "Velocidad máxima")

Que otras cosas podemos mapear?

names(calles)
##  [1] "osm_id"               "name"                 "access"              
##  [4] "addr.interpolation"   "area"                 "bench"               
##  [7] "bicycle"              "bicycle_road"         "bridge"              
## [10] "bus"                  "bus.lanes"            "busway.right"        
## [13] "crossing"             "cycleway"             "cycleway.left"       
## [16] "foot"                 "footway"              "gauge"               
## [19] "hazard"               "highway"              "horse"               
## [22] "incline"              "lanes"                "lanes.backward"      
## [25] "lanes.forward"        "layer"                "lit"                 
## [28] "maxspeed"             "maxspeed.conditional" "motor_vehicle"       
## [31] "motor_vehicle.lanes"  "note"                 "oneway"              
## [34] "oneway.bicycle"       "opposite_track"       "public_transport"    
## [37] "railway"              "segregated"           "service"             
## [40] "shelter"              "short_name"           "sidewalk"            
## [43] "sorting_name"         "source"               "source.maxspeed"     
## [46] "surface"              "tactile_paving"       "traffic_calming"     
## [49] "turn.lanes"           "vehicle"              "geometry"

Podemos graficar tambien calles por cantidad de carriles

ggplot(calles) +
    geom_sf(aes(size = lanes)) +

    labs(title = "Parque Patricios",
         subtitle = "Cantidad de carriles",
         caption = "fuente: OpenStreetMap",
         size = "Cantidad de carriles")

Teniendo en cuenta que se trata de una zona logística de galpones y circulación de tráfico pesado, donde pronto aumentará la circulación de peatones, se podría usar el plano de velocidades máximas y ancho de calles para ver donde sería oportuno colocar semaforos, o reductores de velocidad.

Los dos mapas superpuestos:

ggplot(calles) +
    geom_sf(aes(size = lanes)) +
    geom_sf(aes(color = maxspeed), alpha = 0.5) +

    labs(title = "Parque Patricios",
         subtitle = "Velocidad y Cantidad de carriles",
         caption = "fuente: OpenStreetMap",
         color = "Velocidad Maxima",
         size = "Cantidad de carriles")

Agregamos los semaforos

semaforo <- opq(bbox) %>% 
  add_osm_feature(key = "highway", value = "traffic_signals") %>% 
  osmdata_sf()

A diferencia de las calles, para los semaforos queremos los puntos unicamente

semaforo <- st_intersection(semaforo$osm_points, bbox_poly)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

Graficamos el mapa anterior con la superposición de semaforos

ggplot(calles) +
    geom_sf(aes(size = lanes)) +
    geom_sf(aes(color = maxspeed), alpha = 0.5) +
    geom_sf(data = semaforo, 
            color = "red")+

    labs(title = "Parque Patricios",
         subtitle = "Velocidad y Cantidad de carriles",
         caption = "fuente: OpenStreetMap",
         color = "Velocidad Maxima",
         size = "Cantidad de carriles")

Se nota una concentración de semaforos sobre las avenidas, Un buen lugar para agregar semáforo podría ser en la Av. Amancio Alcorta, que no cuenta con ninguno en el tramo de +400m entre Lavardén y Colonia.

ggplot(calles) +
    geom_sf(aes(size = lanes)) +
    geom_sf(aes(color = maxspeed), alpha = 0.5) +
    geom_sf(data = semaforo, 
            color = "red")+
    geom_sf_label(data = calles%>%
                  filter (name == "Lavardén" | name == "Avenida Colonia"), 
                  aes(label = name), size = 2)+

    labs(title = "Parque Patricios",
         subtitle = "Velocidad y Cantidad de carriles",
         caption = "fuente: OpenStreetMap",
         color = "Velocidad Maxima",
         size = "Cantidad de carriles")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
## not give correct results for longitude/latitude data

Tambien podemos traer informacion de equipamientos y amenidades urbanas, veamos los cafes

cafe <- opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "cafe") %>% 
  osmdata_sf()

A diferencia de las calles, para los Cafes queremos los puntos unicamente

cafe <- st_intersection(cafe$osm_points, bbox_poly)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

Y ahora podemos mapear

ggplot() +
  geom_sf(data = calles, 
            color = "darkslateblue") +
  geom_sf(data = cafe, 
            color = "red")

La oferta es bastante restringida, y se concentra en su mayoría entorno al Parque de los Patricios y la avenida Brasil. Si bien el proyecto de Estacion Buenos Aires contempla 73 locales comerciales, sería bueno que la oferta se atomize en el entorno para que la gente ´salga´ del complejo y evitar el efecto de aghetización que se produjo en otros complejos como Lugano I y II.

La info de OMS tiene otros atributos de interes como el nombre del local, si tienen delivery, si son franquicia, mesas en la vereda y horarios de cierre.

ggplot() +
  geom_sf(data = calles, 
            color = "darkslateblue") +
  geom_sf(data = cafe, 
            aes(color = takeaway)) +
  geom_sf_label(data = cafe, 
                     aes(label = name), size = 3)
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may
## not give correct results for longitude/latitude data

Por último, vamos a buscar un atributo que podamos comparar con la base de datos de BA DATA, para contrastar la información. Vamos a epxlorar la info de estaciones de servicio que se encuentra disponible tanto en OSM como en BA Data.

ES <- opq(bbox) %>% 
  add_osm_feature(key = "amenity", value = "fuel") %>% 
  osmdata_sf()

aislamos

ES <- st_intersection(ES$osm_points, bbox_poly)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

mapeamos

ggplot() +
  geom_sf(data = calles, 
            color = "darkslateblue") +
  geom_sf(data = ES, 
            color = "red")

A pesar de que son pocas, vamos a crear un contador

ES.c <- ES %>% 
    summarise(cantidad = n())
summary(ES.c)
##     cantidad           geometry
##  Min.   :17   MULTIPOINT   :1  
##  1st Qu.:17   epsg:4326    :0  
##  Median :17   +proj=long...:0  
##  Mean   :17                    
##  3rd Qu.:17                    
##  Max.   :17

Por algun motivo, algunas estaciones de servicio se estan contando mas de una vez. Se puede ver tanto en el mapa (las que se suporponen aparecen mas grandes) como en el conteo, donde aparecen 17 en lugar de 6.

Veamos la base de datos de BA

barrios <- st_read('http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios.geojson')
## Reading layer `barrios_badata' from data source `http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios.geojson' using driver `GeoJSON'
## Simple feature collection with 48 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.70529 xmax: -58.33515 ymax: -34.52649
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ES.BA <- read.csv('C:/MEU (DiTella)/Ciencia de datos 2/BA/ES/estaciones_servicio_caba.csv', encoding = "UTF-8")
ggplot() +
  geom_sf(data = barrios) +
  geom_point(data = ES.BA, 
               aes(x = long, y = lat),
               alpha = .3, 
               color = "red")

Filtramos

ES.BA <- ES.BA %>% 
    filter(!is.na(long), !is.na(lat)) %>% 
    st_as_sf(coords = c("long", "lat"), crs = 4326)

Realizamos el join espacial

ES.barrios <- st_join(ES.BA, barrios)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Creamos un contador de Estaciones de Servicio por barrio

ES.barrio.c <- ES.barrios %>% 
    group_by(barrio.x) %>% 
    summarise(cantidad = n())

Visualizamos datos

ggplot(ES.barrio.c) +
    geom_bar(aes(x = barrio.x, weight = cantidad))+
    coord_flip()

Vamos a filtrar los datos para quedarnos unicamente con los de Parque Patricios

seleccion <- filter(ES.barrio.c, barrio.x == "Parque Patricios")
summary(seleccion)
##              barrio.x    cantidad          geometry
##  Parque Patricios:1   Min.   :4   MULTIPOINT   :1  
##                  :0   1st Qu.:4   epsg:4326    :0  
##  Agronomia       :0   Median :4   +proj=long...:0  
##  Almagro         :0   Mean   :4                    
##  Balvanera       :0   3rd Qu.:4                    
##  Barracas        :0   Max.   :4                    
##  (Other)         :0

Parque Patricios cuenta con 4 Estaciones de Servicio

Grafiquemos las estaciones de servicio de ambos dataset en un mismo plano, verde las de BA Datos, y Rojo las de OSM.

ggplot() +
  geom_sf(data = calles, 
            color = "darkslateblue") +
  geom_sf(data = ES, 
            color = "red", size = 5)+
  geom_sf(data = seleccion, 
            color = "green")

Estarian faltando 2 estaciones de servicio en el dataset de BA Datos

Revisamos en google street view para ver si efectivamente existian esas estaciones

GNC Alberti

Dragon Combustibles

¿Por que faltan Estaciones de Servicios en el dataset oficial de BA? La última actualizacion del archivo de BA Data es de mayo del 2019 y las estaciones parecen viejas por lo que debería ser otro motivo.

TP 3: Capturando y explorando datos de Twitter

Vamos a buscar tweets que mencionen alguno de los términos “parque”, “plaza”, “arbolado” o “arboles” en un radio de 10 millas (~16 km) en torno obelisco de la Ciudad de Buenos Aires.

El parámetro n = 5000 es para limitar la búsqueda a los primeros 5000 tweets hallados.

El parámetro include_rts = FALSE descarta retweets

tweets.BA <- search_tweets(q = "parque OR plaza OR arbolado OR arboles",
              geocode = "-34.603722,-58.381592,10mi",
              include_rts = FALSE,
              n = 5000,
              retryonratelimit = TRUE)
## Registered S3 method overwritten by 'openssl':
##   method      from
##   print.bytes Rcpp

Exploramos el contenido

names(tweets.BA)
##  [1] "user_id"                 "status_id"              
##  [3] "created_at"              "screen_name"            
##  [5] "text"                    "source"                 
##  [7] "display_text_width"      "reply_to_status_id"     
##  [9] "reply_to_user_id"        "reply_to_screen_name"   
## [11] "is_quote"                "is_retweet"             
## [13] "favorite_count"          "retweet_count"          
## [15] "quote_count"             "reply_count"            
## [17] "hashtags"                "symbols"                
## [19] "urls_url"                "urls_t.co"              
## [21] "urls_expanded_url"       "media_url"              
## [23] "media_t.co"              "media_expanded_url"     
## [25] "media_type"              "ext_media_url"          
## [27] "ext_media_t.co"          "ext_media_expanded_url" 
## [29] "ext_media_type"          "mentions_user_id"       
## [31] "mentions_screen_name"    "lang"                   
## [33] "quoted_status_id"        "quoted_text"            
## [35] "quoted_created_at"       "quoted_source"          
## [37] "quoted_favorite_count"   "quoted_retweet_count"   
## [39] "quoted_user_id"          "quoted_screen_name"     
## [41] "quoted_name"             "quoted_followers_count" 
## [43] "quoted_friends_count"    "quoted_statuses_count"  
## [45] "quoted_location"         "quoted_description"     
## [47] "quoted_verified"         "retweet_status_id"      
## [49] "retweet_text"            "retweet_created_at"     
## [51] "retweet_source"          "retweet_favorite_count" 
## [53] "retweet_retweet_count"   "retweet_user_id"        
## [55] "retweet_screen_name"     "retweet_name"           
## [57] "retweet_followers_count" "retweet_friends_count"  
## [59] "retweet_statuses_count"  "retweet_location"       
## [61] "retweet_description"     "retweet_verified"       
## [63] "place_url"               "place_name"             
## [65] "place_full_name"         "place_type"             
## [67] "country"                 "country_code"           
## [69] "geo_coords"              "coords_coords"          
## [71] "bbox_coords"             "status_url"             
## [73] "name"                    "location"               
## [75] "description"             "url"                    
## [77] "protected"               "followers_count"        
## [79] "friends_count"           "listed_count"           
## [81] "statuses_count"          "favourites_count"       
## [83] "account_created_at"      "verified"               
## [85] "profile_url"             "profile_expanded_url"   
## [87] "account_lang"            "profile_banner_url"     
## [89] "profile_background_url"  "profile_image_url"

Obtenemos un top 5 de los tweets más populares (con más retweets), su procedencia, y el contenido del tweet:

tweets.BA %>% 
    top_n(5, retweet_count) %>% 
    arrange(desc(retweet_count)) %>% 
    select(screen_name, retweet_count, location, text)
## # A tibble: 5 x 4
##   screen_name   retweet_count location        text                         
##   <chr>                 <int> <chr>           <chr>                        
## 1 MagaliNarvaja          7262 Buenos Aires, ~ unas ganas de tener a una pe~
## 2 Lautiroman19~          3272 Buenos Aires, ~ AHORA. Impresionante. Plaza ~
## 3 NathaRamirez~          2343 Buenos Aires, ~ "Hace unos meses arreglé por~
## 4 WolffWaldo             1935 Buenos Aires, ~ "El  #8SApoyamosALaLeona por~
## 5 zlotomarcelo           1695 Ciudad Autónom~ El macrismo repartiendo merc~

El tweet de magalí buscando gente nueva fue el mas popular. De los 4 puestos siguientes, 3 eran de carácter político.

Por otro lado podemos ver el horario del dia en que se realiza la mayor cantidad de tweets, los diferentes dias de la semana

ts_plot(tweets.BA, "days")

El pico de actividad (dentro de la muestra que tomamos) se da el Sabado 31/8 y el domingo 8/9.

tweets.dia <- tweets.BA %>%
  mutate(dia = as.Date(substr(created_at, 1,10))) %>%
  filter(dia == as.Date("2019-08-10"))

Para ver la popularidad de los usuarios, obtenemos un top 5 de los usuarios más populares (con más seguidores), su procedencia, y el contenido del tweet:

tweets.BA %>% 
    top_n(5, followers_count) %>% 
    arrange(desc(followers_count)) %>% 
    select(screen_name, followers_count, location, text)
## # A tibble: 6 x 4
##   screen_name  followers_count location       text                         
##   <chr>                  <int> <chr>          <chr>                        
## 1 mauriciomac~         4936435 Buenos Aires,~ En vivo desde Córdoba: el Pr~
## 2 clarincom            2911927 Buenos Aires,~ María Eugenia Vidal tendrá s~
## 3 clarincom            2911927 Buenos Aires,~ "Tragedia en el parque nacio~
## 4 clarincom            2911927 Buenos Aires,~ Tras dos años de reclamo vec~
## 5 clarincom            2911927 Buenos Aires,~ Tras cinco meses de obras in~
## 6 clarincom            2911927 Buenos Aires,~ Sigue la polémica en Tierra ~

En este caso, a diferencia del top5 de retweets, los 5 primeros son de caracter político.

Veamos la distribución las cuentas mas populares

 options(scipen = 20)

ggplot(tweets.BA)+
    geom_histogram(aes(x = followers_count))+
     labs(title = "Distribución de popularidad",
       x = "Cantidad followers",
       y = "Cuentas") +
  theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Se ve como la gran mayoria de usuarios tiene pocos o nada de followers, y unos pocos ´Tweetstars´ concentran muchisimos.

Este tipo de distribución se ve en otras redes sociales, como YouTube, donde unos pocos usuarios generan los contenidos consumidos por la mayoria. Recientemente esta minoría, los Youtubers, se unieron al sindicato de comercio mas grande de Europa y exigen a la empresa ser tratados como socios.

FairTube: The YouTubers Union Is Not Messing Around

Vamos a aislar los tweets que poseen coordenadas geográficas (lat y long), crear mapas que muestren posición de los tweets y cantidad de seguidores del usuario que tuitea.

Primero creamos una función que extrae los valores de lat y long del campo problemático

coordenadas <- function(campo_coordenadas) {
    extraer_coordenadas <- function(lista_coords) {
        data_frame(lon = lista_coords[1],
                   lat = lista_coords[2])
    }
    
    map_df(campo_coordenadas, extraer_coordenadas)
}

Y con una cadena de funciones extraemos los datos de georeferenciamiento, los agregamos al dataframe en sus propias columnas, y por último descartamos los campos en formato lista:

tweets_EEVV <- tweets.BA %>% 
    cbind(coordenadas(tweets.BA$coords_coords)) %>% 
    select(-geo_coords, -coords_coords, -bbox_coords) 
## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.

Antes de visualizar, filtramos nuestros tweets para conservar sólo los que contienen coordenadas exactas de posición.

tweets_EEVV_geo <- tweets_EEVV %>% 
    filter(!is.na(lat), !is.na(lon))

El resultado evidencia que los tweets georeferenciados son sólo una pequeña fracción del total que se produce:

nrow(tweets_EEVV_geo)
## [1] 325

Delimitamos los bordes externos del boundinbox

bbox <- c(min(tweets_EEVV_geo$lon),
          min(tweets_EEVV_geo$lat),
          max(tweets_EEVV_geo$lon),
          max(tweets_EEVV_geo$lat))
require(leaflet)
gama <- colorNumeric(palette = "plasma",
                        domain = tweets_EEVV_geo$followers_count)
leaflet(tweets_EEVV_geo) %>% 
    addTiles() %>% 
    addCircleMarkers(popup = ~text,
                     color = ~gama(followers_count)) %>% 
    addLegend(title = "Popularidad", pal = gama, values = ~followers_count)
## Assuming "lon" and "lat" are longitude and latitude, respectively

El tweet del usuario mas popular, un runner que nos regala la siguiente postal desde el Parque 3 de Febrero: "Pasadas en la mañana de Palermo. Los días no pueden ser más lindos, espectacular el cielo del Parque. Mi gente a pleno. Tarea cumplida. #skechersargentina #correrayuda #correrparavivir… https://t.co/ps4Bn7VysJ"

###TP4 - Descubriendo patrones temporales y espaciales en los datos

Vamos a trabajar con los datos de SUACI (sist. unico de atencion ciudadana), para ver los patrones espacio temporales de las quejas y reclamos de los vecinos de la Ciudad de Buenos Aires. Si bien hay datos del 2019, vamos a trabajar con los del año pasado, para poder observar variaciones a lo largo de todo el año calendario.

suaci <- read.csv("C:/MEU (DiTella)/2019 - A2/2 Trimestre/MU115 - Ciencia de datos 2/Datos2/BA/SUACI/sistema-unico-de-atencion-ciudadana-2018.csv", encoding="UTF-8")

Exploramos los datos

dim(suaci)
## [1] 893291     19
names(suaci)
##  [1] "contacto"                  "periodo"                  
##  [3] "categoria"                 "subcategoria"             
##  [5] "concepto"                  "tipo_prestacion"          
##  [7] "fecha_ingreso"             "hora_ingreso"             
##  [9] "domicilio_cgpc"            "domicilio_barrio"         
## [11] "domicilio_calle"           "domicilio_altura"         
## [13] "domicilio_esquina_proxima" "lat"                      
## [15] "long"                      "canal"                    
## [17] "genero"                    "estado_del_contacto"      
## [19] "fecha_cierre_contacto"
head(suaci)
##      contacto periodo                  categoria        subcategoria
## 1 00000001/18  201801                   TRÁNSITO       DENUNCIA VIAL
## 2 00000002/18  201801                   TRÁNSITO       DENUNCIA VIAL
## 3 00000003/18  201801                   TRÁNSITO       DENUNCIA VIAL
## 4 00000004/18  201801 ARBOLADO Y ESPACIOS VERDES PLANTACIÓN DE ÁRBOL
## 5 00000005/18  201801                   TRÁNSITO       DENUNCIA VIAL
## 6 00000006/18  201801                   TRÁNSITO       DENUNCIA VIAL
##                   concepto tipo_prestacion fecha_ingreso  hora_ingreso
## 1 VEHÍCULO MAL ESTACIONADO        DENUNCIA    2018-01-01 12:03:30 a.m.
## 2 VEHÍCULO MAL ESTACIONADO        DENUNCIA    2018-01-01 12:19:33 a.m.
## 3 VEHÍCULO MAL ESTACIONADO        DENUNCIA    2018-01-01 12:19:33 a.m.
## 4      PLANTACIÓN DE ÁRBOL       SOLICITUD    2018-01-01 12:19:36 a.m.
## 5 VEHÍCULO MAL ESTACIONADO        DENUNCIA    2018-01-01 12:20:33 a.m.
## 6 VEHÍCULO MAL ESTACIONADO        DENUNCIA    2018-01-01 12:20:58 a.m.
##   domicilio_cgpc domicilio_barrio       domicilio_calle domicilio_altura
## 1      COMUNA 13            NUÑEZ        PAZ, GRAL. AV.             1020
## 2      COMUNA 13            NUÑEZ                 VEDIA             2225
## 3      COMUNA 13            NUÑEZ    VUELTA DE OBLIGADO             4759
## 4      COMUNA 14          PALERMO REPUBLICA DE LA INDIA             2866
## 5      COMUNA 13            NUÑEZ                 VEDIA             2225
## 6       COMUNA 6        CABALLITO             GUAYAQUIL              325
##   domicilio_esquina_proxima       lat      long canal    genero
## 1                           -34.53893 -58.47435   App masculino
## 2                           -34.53902 -58.47312   App masculino
## 3                           -34.53969 -58.47343   App masculino
## 4                           -34.57929 -58.41552   App masculino
## 5                           -34.53902 -58.47312   App masculino
## 6                           -34.62006 -58.43336   App  femenino
##   estado_del_contacto fecha_cierre_contacto
## 1             Cerrado            2018-01-03
## 2             Cerrado            2018-01-01
## 3             Cerrado            2018-01-03
## 4             Cerrado            2018-09-05
## 5             Cerrado            2018-01-03
## 6             Cerrado            2018-01-03

El dataset cuenta con la informacion necesaria para el ejercicio: fechas y coordenadas.

Ordenamos las categorias con mas reclamos.

suaci %>% 
    count(categoria) %>% 
    top_n(5) %>% 
    arrange(desc(n))
## Selecting by n
## # A tibble: 5 x 2
##   categoria                       n
##   <fct>                       <int>
## 1 LIMPIEZA Y RECOLECCIÓN     296730
## 2 TRÁNSITO                   288098
## 3 ARBOLADO Y ESPACIOS VERDES  69348
## 4 CALLES Y VEREDAS            65094
## 5 ALUMBRADO                   39781

Y analizamos la temporalidad de los contactos relacionados a Arbolado y Espacios Verdes.

suaci %>% 
    filter(categoria == "ARBOLADO Y ESPACIOS VERDES") %>% 
    ggplot() +
        geom_bar(aes(x = month(fecha_ingreso, label = TRUE)))+
  
  labs(title = "Reclamos vinculados a Arbolado y EE.VV",
       caption = "Fuente: BA Data",
       y = "cantidad",
       x= "mes")

A simple vista se puede observar una tendencia alcista hacia el verano

Repetimos el gráfico, pero esta vez incluimos las subcategorias para ver si podemos explicar la variacion mensual de los reclamos.

suaci %>% 
    filter(categoria == "ARBOLADO Y ESPACIOS VERDES") %>% 
    ggplot() +
        geom_bar(aes(x = month(fecha_ingreso, label = TRUE), fill = subcategoria))+
  
  labs(title = "Reclamos vinculados a Arbolado y EE.VV",
       subtitle = "Según tipo de reclamos", 
       caption = "Fuente: BA Data",
       fill = "Subcategorias",
       y = "cantidad",
       x= "mes")

Casi toda la variación temporal se explica en la categoría “Problema con intervención de Arbolado”. Es la unica subcategoría con estacionalidad fuerte.

Para ver que pasa con los demas reclamos vamos a retomar el top 5

suaci_frecuentes <- suaci %>% 
    count(categoria) %>% 
    top_n(5) %>% 
    pull(categoria)
## Selecting by n

Y graficamos en forma de lineas

# Primero realizamos un conteo de delitos por tipo y por mes del año
conteo <-  suaci %>% 
    filter(categoria %in% suaci_frecuentes) %>% 
    count(categoria, mes = month(fecha_ingreso, label = TRUE))
   
# Y ahora a mostras las cantidades mensuales como líneas 
ggplot(conteo) +
    geom_line(aes(x = mes, y = n, group = categoria, color = categoria))+
  expand_limits(y = 0)+
      theme_minimal()

Podemos ver rapidamente que los reclamos por recolección y los de transito van muy por encima del resto de las categorias. Tambien podemos ver que la categoria de transito tiene una estacionalidad mucho mas marcada que las demas: los reclamos por transito bajan mucho los meses de vacaciones, cuando estan cerrados los colegios y la mayoria de las familias se toman vacaciones.

Veamos como se comportan los reclamos durante la semana

# Primero realizamos un conteo de delitos por tipo y por mes del año
conteo <-  suaci %>% 
    filter(categoria %in% suaci_frecuentes) %>% 
    count(categoria, diasemana = wday(fecha_ingreso, label = TRUE))
   
# Y ahora a mostras las cantidades mensuales como líneas 
ggplot(conteo) +
    geom_line(aes(x = diasemana, y = n, group = categoria, color = categoria))+
  expand_limits(y = 0)+
      theme_minimal()

En todas las categorias se observa una variacion importante entre el domingo y el lunes, y entre el viernes y el sabado.

###Analisis Espacial

Ya cargamos la libreria ggmap, ahora vamos a procurarnos un mapa de Buenos Aires. Primero filtramos los datos sin coordenadas, luego obtenemos el “bounding box” de nuestros datos, y finalmente se los pasamos a get_stamenmap().

suaci <- suaci %>% 
    filter(lat <0, long <0)

bbox <- c(min(suaci$long, na.rm = TRUE),
          min(suaci$lat, na.rm = TRUE),
          max(suaci$long, na.rm = TRUE),
          max(suaci$lat, na.rm = TRUE))

CABA <- get_stamenmap(bbox = bbox, 
                      maptype = "toner-lite",
                      zoom=13)
## Source : http://tile.stamen.com/toner-lite/13/2764/4934.png
## Source : http://tile.stamen.com/toner-lite/13/2765/4934.png
## Source : http://tile.stamen.com/toner-lite/13/2766/4934.png
## Source : http://tile.stamen.com/toner-lite/13/2767/4934.png
## Source : http://tile.stamen.com/toner-lite/13/2768/4934.png
## Source : http://tile.stamen.com/toner-lite/13/2764/4935.png
## Source : http://tile.stamen.com/toner-lite/13/2765/4935.png
## Source : http://tile.stamen.com/toner-lite/13/2766/4935.png
## Source : http://tile.stamen.com/toner-lite/13/2767/4935.png
## Source : http://tile.stamen.com/toner-lite/13/2768/4935.png
## Source : http://tile.stamen.com/toner-lite/13/2764/4936.png
## Source : http://tile.stamen.com/toner-lite/13/2765/4936.png
## Source : http://tile.stamen.com/toner-lite/13/2766/4936.png
## Source : http://tile.stamen.com/toner-lite/13/2767/4936.png
## Source : http://tile.stamen.com/toner-lite/13/2768/4936.png
## Source : http://tile.stamen.com/toner-lite/13/2764/4937.png
## Source : http://tile.stamen.com/toner-lite/13/2765/4937.png
## Source : http://tile.stamen.com/toner-lite/13/2766/4937.png
## Source : http://tile.stamen.com/toner-lite/13/2767/4937.png
## Source : http://tile.stamen.com/toner-lite/13/2768/4937.png
## Source : http://tile.stamen.com/toner-lite/13/2764/4938.png
## Source : http://tile.stamen.com/toner-lite/13/2765/4938.png
## Source : http://tile.stamen.com/toner-lite/13/2766/4938.png
## Source : http://tile.stamen.com/toner-lite/13/2767/4938.png
## Source : http://tile.stamen.com/toner-lite/13/2768/4938.png

Para empezar vamos a mapear densidades que muestre donde se concentran la mayor cantidad de observaciones.

ggmap(CABA) +
    geom_bin2d(data = suaci, aes(x = long, y = lat), bins = 75) +
    scale_fill_viridis_c()
## Warning: Removed 4 rows containing missing values (geom_tile).

Se ve que la mayor cantidad de reclamos se centra en el corredor Norte, pero no es muy claro el mapa. Vamos a probar con un kernel density

ggmap(CABA) +
    geom_density2d(data = suaci, aes(x = long, y = lat, color = stat(level))) +
    scale_color_viridis_c(option = "D", direction = -1)

Podemos ver que hay varios puntos calientes a lo largo del corredor norte, y un foco de reclamos en Barracas, cerca del acceso a Capital por Puente Pueyrredon (congestión vehicular?)

Si empezamos a visualizar por categorias podemos darnos una idea mas acabada de que pasa en esos hotspots. Tomamos las 5 categorias del Top5.

ggmap(CABA) +
    geom_point(data = filter(suaci, categoria %in% suaci_frecuentes),
                      aes(x = long, y = lat, color = categoria),
                      size = 0.02, alpha = 0.1)+
   guides(color = guide_legend(override.aes = list(size=2, alpha = 1))) +
    scale_color_brewer(palette = "Set1")

Por la cantidad de puntos, la info resulta dificil de leer, asi que armamos un mapa de facetado por categorias

ggmap(CABA) +
    geom_point(data = filter(suaci, categoria %in% suaci_frecuentes), 
               aes(x = long, y = lat, color = categoria),
               size = 0.005, alpha = 0.1) +
    scale_color_brewer(palette = "Set1") +
    facet_wrap(~categoria)+
  
  guides(color = guide_legend(override.aes = list(size=2, alpha = 1))) +
    scale_color_brewer(palette = "Set1")+
 
    labs(title = "Distribución geográfica por tipo de reclamos",
         color = "Top 5 reclamos")
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Para hacer las diferencias aún mas nítidas, podemos facetar una estimación de densidad:

ggmap(CABA) +
    geom_density2d(data = filter(suaci, categoria %in% suaci_frecuentes), 
                   aes(x = long, y = lat, color = stat(level))) +
    scale_color_viridis_c(option = "D", direction = -1) +
    facet_wrap(~categoria)+
  
      labs(title = "Densidad geográfica por tipo de reclamos",
         color = "Top 5 reclamos")

En el caso del tránsito se observa una concentración de reclamos en la zona norte, mientras que en limpieza el foco se traslada a Recoleta. Arbolado es la mas homogénea en terminos de localización de reclamos.

Veamos que pasa geográficamente con los reclamos de limpieza y recolección a lo largo del año:

suaci <- suaci %>% 
    mutate(mes = month(fecha_ingreso, label = TRUE))


ggmap(CABA) +
    geom_density2d(data = filter(suaci, 
                                 categoria == "LIMPIEZA Y RECOLECCIÓN"),
                   aes(x = long, 
                       y = lat, 
                       color = stat(level))) +
  
    scale_color_viridis_c(option = "D", direction = -1) +
    facet_wrap(~mes, nrow = 3) +
    labs(title = "Concentración espacial de reclamos de limpieza",
         subtitle = "según el mes")

Se puede ver como los patrones espaciales mantienen ciertas similaridades a lo largo de los meses (concentración en la zona de recoleta), siendo mayo uno de los meses más atípicos, con un hotspot de reclamos en la zona sur.

TP 5: Analizando y visualizando flujos de viajes urbanos

Vamos a analizar los viajes origen-destino del sistema de EcoBici 2018

bici <- read.csv('C:/MEU (DiTella)/2019 - A2/2 Trimestre/MU115 - Ciencia de datos 2/Datos2/BA/bici/recorridos-realizados-2018.csv', encoding = "UTF-8")

Exploramos el dataset

dim(bici)
## [1] 2619968       9
names(bici)
## [1] "bici_id_usuario"              "bici_Fecha_hora_retiro"      
## [3] "bici_tiempo_uso"              "bici_nombre_estacion_origen" 
## [5] "bici_estacion_origen"         "bici_nombre_estacion_destino"
## [7] "bici_estacion_destino"        "bici_sexo"                   
## [9] "bici_edad"
summary(bici)
head(bici)
##   bici_id_usuario bici_Fecha_hora_retiro           bici_tiempo_uso
## 1            5453    2018-01-01 00:08:05 0 days 00:19:53.000000000
## 2             673    2018-01-01 00:18:05 0 days 00:26:19.000000000
## 3          179119    2018-01-01 00:20:14 0 days 00:27:39.000000000
## 4          400147    2018-01-01 00:20:22 0 days 00:48:51.000000000
## 5          400156    2018-01-01 00:20:31 0 days 00:49:27.000000000
## 6          476733    2018-01-01 00:21:01 0 days 00:36:10.000000000
##   bici_nombre_estacion_origen bici_estacion_origen
## 1                     Uruguay                   45
## 2                     Posadas                  189
## 3          Hospital Rivadavia                   50
## 4              Macacha Güemes                  111
## 5              Macacha Güemes                  111
## 6                       Yatay                  121
##    bici_nombre_estacion_destino bici_estacion_destino bici_sexo bici_edad
## 1               Virrey Cevallos                   183         M        45
## 2                 Guardia Vieja                   110         M        61
## 3                       Padilla                    31         F        52
## 4             Acuña de Figueroa                    54         M        27
## 5             Acuña de Figueroa                    54         F        27
## 6 Billinghurst y Valentin Gomez                   143         F        31

Tenemos la cantidad de viajes, pero queremos conocer la cantidad de estaciones, por lo que creamos un contador

conteo.est <- bici %>% 
    group_by(bici_nombre_estacion_origen) %>% 
    summarise(cantidad = n())
dim(conteo.est)
## [1] 199   2

Ya tenemos los recorridos y cantidad de estaciones, ahora vamos a descargar el dataset de Estaciones para ver la ubicacion, que es el dato que nos esta faltando.

estaciones <- read.csv('C:/MEU (DiTella)/2019 - A2/2 Trimestre/MU115 - Ciencia de datos 2/Datos2/BA/bici/estaciones-de-bicicletas-publicas.csv', encoding = "UTF-8")

Exploramos

names(estaciones)
##  [1] "long"      "lat"       "nombre"    "domicilio" "imagen"   
##  [6] "automat"   "observ"    "nro_est"   "horario"   "dire_norm"
dim(estaciones)
## [1] 199  10

Ambos datasets nos indican que hay 199 estaciones. Vamos a visualizarlas

Definimos el BoundingBox

library(ggmap)

bbox <- make_bbox(estaciones$X, estaciones$Y)
## Warning in min(x, na.rm = na.rm): ningún argumento finito para min;
## retornando Inf
## Warning in max(x, na.rm = na.rm): ningun argumento finito para max;
## retornando -Inf
## Warning in min(x, na.rm = na.rm): ningún argumento finito para min;
## retornando Inf
## Warning in max(x, na.rm = na.rm): ningun argumento finito para max;
## retornando -Inf
bbox
##   left bottom  right    top 
##   -Inf   -Inf    Inf    Inf
bbox <- c(min(estaciones$lon),
          min(estaciones$lat),
          max(estaciones$lon),
          max(estaciones$lat))
mapa_base <- get_stamenmap(bbox, color = "bw", zoom = 13)
## Source : http://tile.stamen.com/terrain/13/2765/4935.png
## Source : http://tile.stamen.com/terrain/13/2766/4935.png
## Source : http://tile.stamen.com/terrain/13/2767/4935.png
## Source : http://tile.stamen.com/terrain/13/2768/4935.png
## Source : http://tile.stamen.com/terrain/13/2765/4936.png
## Source : http://tile.stamen.com/terrain/13/2766/4936.png
## Source : http://tile.stamen.com/terrain/13/2767/4936.png
## Source : http://tile.stamen.com/terrain/13/2768/4936.png
## Source : http://tile.stamen.com/terrain/13/2765/4937.png
## Source : http://tile.stamen.com/terrain/13/2766/4937.png
## Source : http://tile.stamen.com/terrain/13/2767/4937.png
## Source : http://tile.stamen.com/terrain/13/2768/4937.png
ggmap(mapa_base) +
     geom_point(data = estaciones, 
               aes(x = long, y = lat),
               alpha = .8, 
               color = "orange")

Si queremos saber cuales son los trayectos mas populares, podemos armar un contador que contabilize pares de estaciones

conteo <- bici %>% 
    group_by(bici_estacion_origen, bici_estacion_destino) %>% 
    summarise(total = n())

Vamos a hacer un heatmap para evaluar el grado de interconexion, que muestre la cantidad de viaje entre pares.

ggplot() + 
    geom_tile(data = conteo, aes(x = bici_estacion_origen, y = bici_estacion_destino, fill = total)) +
    scale_fill_distiller(palette = "Spectral")

La numeración discontinua de la estaciones dificulta la lectura del mismo.

Para evitar el bache en el mapa de calor vamos a tratar a las estaciones como una variable categórica (un factor) en lugar de numérica.

ggplot() + 
    geom_tile(data = conteo, 
              aes(x = as.factor(bici_estacion_origen),
                  y = as.factor(bici_estacion_destino),
                  fill = total)) +
    scale_fill_distiller(palette = "Spectral")

La densidad del grafico hace dificil sacar conclusiones, pero podemos ver en la linea verde que la mayoria de los viajes empiezan y terminan en la misma estación. Podria ser por una cuestion recreativa, el usuario da una vuelta y la deja en el mismo lugar. Otra explicación podria ser que al sacarla ven que no esta en condiciones óptimas y deciden dejarla para agarrar otra. Esto último se podría verificar por ejemplo filtrando los viajes menores a 1 minuto.

Otra información que nos proporciona el gráfico es la concentración de puntos rojos (mucha cantidad) hacia el cero del eje de coordenadas, es decir proximo a la estación cero. Esto indica que estos viajes son en la zona centrica, donde se pusieron las primeras estaciones y por lo tanto la que tienen numero de serie menor.

Vamos a obtener los 30 viajes mas populares, descartando los viajes “circulares” (con el mismo origen y destino):

top30 <- conteo %>% 
    ungroup() %>% 
    filter(bici_estacion_origen != bici_estacion_destino) %>% 
    top_n(30)
## Selecting by total
top30
## # A tibble: 30 x 3
##    bici_estacion_origen bici_estacion_destino total
##                   <int>                 <int> <int>
##  1                    1                     2  1111
##  2                    1                   103  1219
##  3                    1                   131  1408
##  4                    2                     1   932
##  5                    5                    14  1023
##  6                    5                    44  1307
##  7                    5                   160  1151
##  8                    5                   177  2213
##  9                    9                    30   967
## 10                    9                    66  1091
## # ... with 20 more rows
ggplot() + 
    geom_tile(data = top30, 
              aes(x = as.factor(bici_estacion_origen),
                  y = as.factor(bici_estacion_destino),
                  fill = total)) +
    scale_fill_distiller(palette = "Spectral")

y el top 10

top10 <- conteo %>% 
    ungroup() %>% 
    filter(bici_estacion_origen != bici_estacion_destino) %>% 
    arrange(desc(total)) %>% 
    top_n(10)
## Selecting by total
top10
## # A tibble: 10 x 3
##    bici_estacion_origen bici_estacion_destino total
##                   <int>                 <int> <int>
##  1                  177                     5  2323
##  2                    5                   177  2213
##  3                   14                   160  1560
##  4                    1                   131  1408
##  5                  131                     1  1407
##  6                   26                   131  1362
##  7                    5                    44  1307
##  8                  103                     1  1270
##  9                  160                    14  1253
## 10                   14                   159  1232
ggplot() + 
    geom_tile(data = top10, 
              aes(x = as.factor(bici_estacion_origen),
                  y = as.factor(bici_estacion_destino),
                  fill = total)) +
    scale_fill_distiller(palette = "Spectral")

y el top 5

top5 <- conteo %>% 
    ungroup() %>% 
    filter(bici_estacion_origen != bici_estacion_destino) %>% 
    arrange(desc(total)) %>% 
    top_n(5)
## Selecting by total
top5
## # A tibble: 5 x 3
##   bici_estacion_origen bici_estacion_destino total
##                  <int>                 <int> <int>
## 1                  177                     5  2323
## 2                    5                   177  2213
## 3                   14                   160  1560
## 4                    1                   131  1408
## 5                  131                     1  1407

Ahora vamos a mapear estos 5 viajes mas frecuentes

Primero hacemos un join del dataframe con el conteo de viajes contra el de posición de estaciones, para agregar las coordenadas

top10 <- top10 %>% 
    left_join(estaciones[c("long", "lat", "nombre", "nro_est")], 
              by = c("bici_estacion_origen" = "nro_est")) %>% 
    rename(ORIGEN_X = long,
           ORIGEN_Y = lat,
           ORIGEN_NOMBRE = nombre)
top10
## # A tibble: 10 x 6
##    bici_estacion_or~ bici_estacion_d~ total ORIGEN_X ORIGEN_Y ORIGEN_NOMBRE
##                <int>            <int> <int>    <dbl>    <dbl> <fct>        
##  1               177                5  2323    -58.4    -34.6 Planetario   
##  2                 5              177  2213    -58.4    -34.6 Plaza Italia 
##  3                14              160  1560    -58.4    -34.6 Pacífico     
##  4                 1              131  1408    -58.4    -34.6 Facultad de ~
##  5               131                1  1407    -58.4    -34.6 Retiro III   
##  6                26              131  1362    -58.4    -34.6 Juana Manso  
##  7                 5               44  1307    -58.4    -34.6 Plaza Italia 
##  8               103                1  1270    -58.4    -34.6 Malba        
##  9               160               14  1253    -58.4    -34.6 Godoy Cruz y~
## 10                14              159  1232    -58.4    -34.6 Pacífico

y hacemos lo propio con las estaciones de destino

top10 <- top10 %>% 
    left_join(estaciones[c("long", "lat", "nombre", "nro_est")], 
              by = c("bici_estacion_destino" = "nro_est")) %>% 
    rename(DESTINO_X = long,
           DESTINO_Y = lat,
           DESTINO_NOMBRE = nombre)

top10
## # A tibble: 10 x 9
##    bici_estacion_o~ bici_estacion_d~ total ORIGEN_X ORIGEN_Y ORIGEN_NOMBRE
##               <int>            <int> <int>    <dbl>    <dbl> <fct>        
##  1              177                5  2323    -58.4    -34.6 Planetario   
##  2                5              177  2213    -58.4    -34.6 Plaza Italia 
##  3               14              160  1560    -58.4    -34.6 Pacífico     
##  4                1              131  1408    -58.4    -34.6 Facultad de ~
##  5              131                1  1407    -58.4    -34.6 Retiro III   
##  6               26              131  1362    -58.4    -34.6 Juana Manso  
##  7                5               44  1307    -58.4    -34.6 Plaza Italia 
##  8              103                1  1270    -58.4    -34.6 Malba        
##  9              160               14  1253    -58.4    -34.6 Godoy Cruz y~
## 10               14              159  1232    -58.4    -34.6 Pacífico     
## # ... with 3 more variables: DESTINO_X <dbl>, DESTINO_Y <dbl>,
## #   DESTINO_NOMBRE <fct>

Vamos a mapear el recorrido del viaje mas popular: Planetario - Pza. Italia