Hemos llegado al final del módulo, y como cierre nos enfocaremos en una herramienta muy útil a la hora de trabajar con datos geográficos de Ciudades:
Cáclulo de ruteos (distancia y duración) entre 2 o más puntos de interés, que lo haremos a partir del uso del paquete osrm.
Comencemos activando (e instalando en caso de ser necesario) todas las librerías que utilizaremos hoy:
#install.packages("tidyverse")
library(tidyverse)
#install.packages("sf")
library(sf)
#install.packages("ggmap")
library(ggmap)
#install.packages("leaflet")
library(leaflet)
#install.packages("osrm")
library(osrm)
La duración de los recorridos que se realizan a diario en la Ciudad y las distancias entre los diferentes puntos de interés son un tema muy relevante a la hora de analizar las dinámicas urbanas. Por suerte desde R podemos conectarnos al servicio de ruteo de OSRM, un paquete de uso libre, basado en datos de OpenStreetMap, que es muy útil a la hora de calcular distancias (km) y tiempos de viaje (minutos) entre 2 o más puntos georreferenciados. Pueden consultar la documentación de la librería aquí.
El paquete se compone de 4 funciones que son:
En la clase de hoy trabajaremos con la función
osrmRoute()
para analizar los recorridos más
realizados con bicicletas públicas (llamadas Blue Bikes) durante la
primera semana de Mayo 2021 en Boston, EEUU.
Estos datos son públicos y pueden encontrarlos en la Página Oficial de Blue Bikes. Pero en este caso, usaremos 2 datasets (en formato .csv) previamente filtrados para facilitar la manipulación de la información durante la clase. Pueden descargarlos desde los siguientes links:
Recomendación: Al descargarlo, moverlo de la carpeta “Descargas” a la carpeta llamada “data” dentro de la carpeta del Proyecto donde estén trabajando.
Comencemos cargando el dataset de estaciones:
estaciones <- read.csv("data/boston_station_xy.csv",
stringsAsFactors = TRUE)
Y exploremos el contenido:
summary(estaciones)
## station_id station_name station_x
## Min. : 1.0 160 Arsenal : 1 Min. :-71.23
## 1st Qu.:103.8 175 N Harvard St : 1 1st Qu.:-71.12
## Median :214.5 18 Dorrance Warehouse : 1 Median :-71.09
## Mean :246.6 191 Beacon St : 1 Mean :-71.09
## 3rd Qu.:394.2 2 Hummingbird Lane at Olmsted Green: 1 3rd Qu.:-71.06
## Max. :510.0 30 Dane St : 1 Max. :-71.01
## (Other) :362
## station_y
## Min. :42.27
## 1st Qu.:42.34
## Median :42.36
## Mean :42.35
## 3rd Qu.:42.37
## Max. :42.42
##
Las columnas son:
station_id: contiene un ID de estación
station_name: contiene el nombre de la estación
station_x: contiene la longitud (coordenada x) de la estación
station_y: contiene la latitud (coordenada y) de la estación
Aprovechemos que tenemos las coordenadas y hagamos una bounding box con ggmap así podemos descargar un mapa de fondo:
bbox_boston <- make_bbox(estaciones$station_x, estaciones$station_y)
mapa_boston <- get_stamenmap(bbox_boston,
zoom = 12)
Y veamos el resultado:
ggmap(mapa_boston)
Ahora ubiquemos las estaciones sobre el mapa:
ggmap(mapa_boston)+
geom_point(data=estaciones, aes(x=station_x, y=station_y), inherit.aes = FALSE)+
labs(title="Estaciones de Bicicleta",
subtitle="Blue Bike - Boston",
caption="Fuente: https://www.bluebikes.com/system-data")+
theme_void()
Podemos observar que la red de estaciones de bicicletas de Blue Bike se distribuye por “todo Boston”, pero la mayor cantidad se concentra en la zona céntrica.
Ahora carguemos el dataset de viajes en bicicleta de la primera semana de mayo 2021 (1 al 7/05):
recorridos <- read.csv("data/boston_bikes_may21.csv",
stringsAsFactors = TRUE)
dim(recorridos)
## [1] 45673 5
Hay 45.673 registros y 5 columnas. Es decir, que durante la primera semana de mayo hubo 45.673 viajes en bicicletas del sistema Blue Bikes. Investiguemos un poco más:
summary(recorridos)
## start_time start_station_id
## 2021-05-01 18:29:56: 5 Min. : 1.0
## 2021-05-01 21:52:08: 5 1st Qu.: 58.0
## 2021-05-01 15:47:15: 4 Median :108.0
## 2021-05-01 16:24:52: 4 Mean :168.9
## 2021-05-01 17:28:10: 4 3rd Qu.:296.0
## 2021-05-01 18:41:38: 4 Max. :510.0
## (Other) :45647
## start_station_name
## MIT at Mass Ave / Amherst St : 1038
## Central Square at Mass Ave / Essex St : 770
## Charles Circle - Charles St at Cambridge St : 733
## Christian Science Plaza - Massachusetts Ave at Westland Ave: 649
## Beacon St at Massachusetts Ave : 584
## Cross St at Hanover St : 568
## (Other) :41331
## end_station_id
## Min. : 1.0
## 1st Qu.: 56.0
## Median :107.0
## Mean :168.2
## 3rd Qu.:282.0
## Max. :510.0
##
## end_station_name
## MIT at Mass Ave / Amherst St : 1067
## Central Square at Mass Ave / Essex St : 783
## Charles Circle - Charles St at Cambridge St : 748
## Christian Science Plaza - Massachusetts Ave at Westland Ave: 649
## Beacon St at Massachusetts Ave : 585
## Cross St at Hanover St : 564
## (Other) :41277
Podemos ver que la estación en donde más viajes comenzaron y terminaron fue “MIT at Mass Ave / Amherst St”.
Analicemos cual es el recorrido que más se realizó, es decir la agrupación origen-destino que aparece con mayor frecuencia en la base de datos:
viajes <- recorridos %>%
group_by(start_station_id, start_station_name, end_station_id, end_station_name) %>%
summarise(cant_viajes = n())
viajes %>%
arrange(desc(cant_viajes)) %>%
head()
## # A tibble: 6 x 5
## # Groups: start_station_id, start_station_name, end_station_id [6]
## start_station_id start_station_na~ end_station_id end_station_name cant_viajes
## <int> <fct> <int> <fct> <int>
## 1 67 MIT at Mass Ave ~ 179 MIT Vassar St 67
## 2 68 Central Square a~ 178 MIT Pacific St ~ 57
## 3 46 Christian Scienc~ 46 Christian Scien~ 56
## 4 179 MIT Vassar St 67 MIT at Mass Ave~ 56
## 5 67 MIT at Mass Ave ~ 67 MIT at Mass Ave~ 52
## 6 68 Central Square a~ 67 MIT at Mass Ave~ 52
El viaje que más se realizó fue entre la estación 67 (MIT at Mass Ave / Amherst St) y la estación 179 (MIT Vassar St).
Observemos la distribución de la variable “cant_viajes” en un histograma:
ggplot(viajes)+
geom_histogram(aes(x=cant_viajes))
Podemos observar que hay muchos recorridos (origen-destino) que se hicieron pocas veces (menos de 10) y pocos recorridos que se hicieron muchas veces (más de 10).
Veamos los resultados en un gráfico de matriz:
ggplot() +
geom_tile(data = viajes, aes(x = start_station_id, y = end_station_id, fill = cant_viajes))+
scale_fill_distiller(palette = "RdYlGn")
La visualización es bastante difícil de interpretar ya que la mayoría de combinaciones origen-destino tienen menos de 20 viajes, pero se ven algunos detalles como por ejemplo una línea recta a 45°. Esto significa que hay varios viajes circulares, es decir que tienen el origen y destino en la misma estación. Esto puede deberse a viajes recreativos por ejemplo.
Veamos el top 10 de recorridos más realizados (sin tener en cuenta los viajes circulares ya que no vamos a poder rutearlos):
top_10 <- viajes %>%
filter(start_station_id != end_station_id) %>%
arrange(desc(cant_viajes)) %>%
head(10)
top_10
## # A tibble: 10 x 5
## # Groups: start_station_id, start_station_name, end_station_id [10]
## start_station_id start_station_n~ end_station_id end_station_name cant_viajes
## <int> <fct> <int> <fct> <int>
## 1 67 MIT at Mass Ave~ 179 MIT Vassar St 67
## 2 68 Central Square ~ 178 MIT Pacific St ~ 57
## 3 179 MIT Vassar St 67 MIT at Mass Ave~ 56
## 4 68 Central Square ~ 67 MIT at Mass Ave~ 52
## 5 67 MIT at Mass Ave~ 53 Beacon St at Ma~ 48
## 6 67 MIT at Mass Ave~ 68 Central Square ~ 48
## 7 157 Seaport Blvd at~ 47 Cross St at Han~ 47
## 8 9 Commonwealth Av~ 446 700 Commonwealt~ 43
## 9 178 MIT Pacific St ~ 80 MIT Stata Cente~ 42
## 10 471 MIT Carleton St~ 67 MIT at Mass Ave~ 42
Tal como vimos anteriormente, el recorrido que más veces se realizó fue desde la estación “MIT at Mass Ave / Amherst St” hasta la estación “MIT Vassar St” y se registró 67 veces.
Llegó el momento de hacer nuestro primer ruteo. Comencemos con un ejemplo de “ruteo simple”, es decir calcular la distancia y el tiempo entre un punto y otro punto.
Para lograr esto utilizaremos la función osrmRoute()
que
necesita que le asignemos toda la información relacionada al origen y al
destino, es decir que debemos tener como mínimo 6 variables que
indiquen:
Por suerte, tenemos toda esa información pero aún está dividida en 2 datasets y deberíamos unificarla para que sea más simple manipularla. Para esto, agreguemos al “top_10” las coordenadas de origen y las de destino:
top_10 <- top_10 %>%
left_join(estaciones, by=c("start_station_id"="station_id",
"start_station_name"="station_name"))
top_10 <- top_10 %>%
rename(start_x=station_x,
start_y=station_y)
top_10 <- top_10 %>%
left_join(estaciones, by=c("end_station_id"="station_id",
"end_station_name"="station_name"))
top_10 <- top_10 %>%
rename(end_x=station_x,
end_y=station_y)
head(top_10)
## # A tibble: 6 x 9
## # Groups: start_station_id, start_station_name, end_station_id [6]
## start_station_id start_station_na~ end_station_id end_station_name cant_viajes
## <int> <fct> <int> <fct> <int>
## 1 67 MIT at Mass Ave ~ 179 MIT Vassar St 67
## 2 68 Central Square a~ 178 MIT Pacific St ~ 57
## 3 179 MIT Vassar St 67 MIT at Mass Ave~ 56
## 4 68 Central Square a~ 67 MIT at Mass Ave~ 52
## 5 67 MIT at Mass Ave ~ 53 Beacon St at Ma~ 48
## 6 67 MIT at Mass Ave ~ 68 Central Square ~ 48
## # ... with 4 more variables: start_x <dbl>, start_y <dbl>, end_x <dbl>,
## # end_y <dbl>
Ahora si, comencemos con un primer ejemplo de un ruteo simple, es decir de un punto a otro punto. Para esto, quedémonos solo con el recorrido más realizado (luego trabajaremos con el resto):
viaje1 <- top_10 %>%
arrange(desc(cant_viajes)) %>%
head(1)
viaje1
## # A tibble: 1 x 9
## # Groups: start_station_id, start_station_name, end_station_id [1]
## start_station_id start_station_na~ end_station_id end_station_name cant_viajes
## <int> <fct> <int> <fct> <int>
## 1 67 MIT at Mass Ave ~ 179 MIT Vassar St 67
## # ... with 4 more variables: start_x <dbl>, start_y <dbl>, end_x <dbl>,
## # end_y <dbl>
Hagamos un ruteo con osrmRoute()
para poder calcular la
distancia y el tiempo de viaje del recorrido más realizado:
ruteo1 <- osrmRoute(src = c(viaje1$start_station_id, viaje1$start_x, viaje1$start_y),
dst = c(viaje1$end_station_id, viaje1$end_x, viaje1$end_y),
returnclass = "sf",
overview = "full",
osrm.profile = "bike")
Como verán tuvimos que asignar 5 parámetros a la función:
En src tenemos que indicar el origen del ruteo con un identificador del sitio (start_station_id) y sus coordenadas.
En dst tenemos que indicar el destino del ruteo con un identificador del sitio (end_station_id) y sus coordenadas.
En returnclass tenemos que indicar el tipo de objeto que queremos que nos devuelva la función, en este caso, un objeto espacial sf (línea de recorrido).
En overview debemos indicar la calidad con la que se genera la geometría: full para obtener la mayor precisión, simplified para menor precisión y FALSE para obtener solo tiempo y distancia (sin geometría).
En osrm.profile indicaremos el modo de transporte en el que queremos realizar el ruteo: car, bike o foot.
Veamos el resultado que obtuvimos:
ruteo1
## Simple feature collection with 1 feature and 4 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -71.10395 ymin: 42.35483 xmax: -71.09321 ymax: 42.35814
## Geodetic CRS: WGS 84
## src dst duration distance geometry
## 67_179 67 179 5.366667 1.0905 LINESTRING (-71.09321 42.35...
El recorrido más realizado en bicicleta fue desde la estación id 67 a la estación id 179, y se detectó una distancia de 1,09 km con una duración de 5,36 min.
Reflejemos este recorrido en el mapa:
ggmap(mapa_boston)+
geom_point(data=estaciones, aes(x=station_x, y=station_y), inherit.aes = FALSE)+
geom_sf(data=ruteo1, color="red", size=1.5, inherit.aes = FALSE)+
labs(title="Recorrido más Realizado en Bicicleta",
subtitle="Blue Bike - Boston",
caption="Fuente: https://www.bluebikes.com/system-data")+
theme_void()
Se ve que el recorrido es dentro de Cambridge pero hagamos un mapa interactivo que nos permita hacer “zoom” para obtener mayor detalle:
leaflet(ruteo1) %>%
addTiles() %>%
addPolylines()
Probemos ajustar un poco la estética y agregar los puntos de origen y
destino con addCircleMarkers()
:
leaflet(ruteo1) %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addPolylines(color = "red",
label = paste("Distancia:", ruteo1$distance, "|",
"Duración:", ruteo1$duration)) %>%
addCircleMarkers(data=viaje1, ~start_x, ~start_y,
popup = ~start_station_name,
color="white") %>%
addCircleMarkers(data=viaje1, ~end_x, ~end_y,
popup = ~end_station_name,
color="white")
Hasta acá ya logramos tener nuestro primer ruteo realizado. Pero, ¿Qué ocurre si en vez de hacerlo con un único recorrido, necesitamos hacerlo para los 10 más realizados?
La función osrmRoute()
fue diseñada para que solamente
se aplique sobre un registro, es decir que no existe un parámetro que
nos permita ejecutarla sobre un dataset de 2 o más registros. Para
lograr esto con lo que sabemos hasta ahora podríamos georreferenciar
todos los registros por separado y luego unirlos con
rbind()
. Sin embargo, existen diferentes métodos que son
muy útiles si queremos evitar tantos pasos: hoy optaremos por una
función de “mapeo” llamada pmap()
, que se compone por un
vector (en este caso las coordenadas de origen y destino) y por una
función (que tendremos que diseñar a partir de
osrmRoute()
).
Para obtener los resultados deseados, vamos a tener que crear lo que en R le llamamos función y luce así:
nombre <- function(argumentos) { operaciones }
De lado izquierdo de la “flechita” irá el nombre que le asignamos a nuestra función y que luego nos servirá para ejecutarla, al igual que como hacemos con cualquier otra función de R que ya conozcamos.
Luego utilizamos la función llamada function()
donde
entre los paréntesis asignaremos los argumentos o variables que la
función va a utilizar al momento de realizar las operaciones. Se
escriben entre “()” y se separan por “,” si son 2 o más.
Por último, asignamos dentro de los {} todas las operaciones que queremos que se realicen cuando la función sea ejecutada.
Para que se entienda mejor como funciona, apliquemos todo esto sobre el dataset “top_10”.
Comencemos creando nuestra primera función, a la que llamaremos ruteo_funcion:
ruteo_funcion <- function(nombre_origen, lon_origen, lat_origen,
nombre_destino, lon_destino, lat_destino)
{
osrmRoute(src = c(nombre_origen, lon_origen, lat_origen),
dst = c(nombre_destino, lon_destino, lat_destino),
returnclass = "sf",
overview = "full",
osrm.profile = "bike")
}
¿Qué hicimos? ¿Qué usamos como argumento y qué usamos como operación?
El argumento elegido se compuso de los 6
elementos sobre los cuales trabajará nuestra función: nombre de origen,
longitud de origen, latitud de origen, nombre de destino, longitud de
destino, latitud de destino. Nótese que el argumento debe aparecer tanto
entre los () de function()
como dentro de
osrmRoute()
.
La operación asignada fue un ruteo con
osrmRoute()
sobre todos los valores del argumento y se
ajustaron 3 parámetros: returnclass, overview y
osrm.profile.
Tal como se mencionó más arriba, un punto importante a destacar es
que la función pmap()
necesita si o si que le asignemos un
vector para que lo tome como input y repita sobre cada elemento
la función elegida, que en este caso será
ruteo_funcion. Por lo tanto es necesario que
transformemos las 6 columnas elegidas del dataset top_10 a un vector con
la función list()
:
ruteo_top10 <- list(top_10$start_station_name, top_10$start_x, top_10$start_y,
top_10$end_station_name, top_10$end_x, top_10$end_y)
Ahora si, ya estamos en condiciones de utilizar pmap()
para rutear todos los recorridos:
ruteo_top10 <- pmap(ruteo_top10, ruteo_funcion)
Veamos que resultado obtuvimos:
summary(ruteo_top10)
## Length Class Mode
## [1,] 5 sf list
## [2,] 5 sf list
## [3,] 5 sf list
## [4,] 5 sf list
## [5,] 5 sf list
## [6,] 5 sf list
## [7,] 5 sf list
## [8,] 5 sf list
## [9,] 5 sf list
## [10,] 5 sf list
Podemos observar que la función nos devolvió una lista con todos los
objetos sf resultantes de los 10 ruteos. Para poder analizar y
mapear los resultados, necesitamos utilizar reduce()
con el
fin de reducir la lista e ir uniendo todas las filas para que nos
devuelva el resultado como un dataset único:
ruteo_top10 <- ruteo_top10 %>%
reduce(rbind)
Veamos que ocurrió:
summary(ruteo_top10)
## src dst duration distance
## Min. : 84.0 Min. : 11.0 Min. :4.028 Min. :0.8867
## 1st Qu.:138.2 1st Qu.: 93.5 1st Qu.:5.110 1st Qu.:0.9417
## Median :232.0 Median :231.0 Median :5.431 Median :1.1229
## Mean :196.1 Mean :165.1 Mean :5.632 Mean :1.1059
## 3rd Qu.:234.5 3rd Qu.:233.2 3rd Qu.:6.119 3rd Qu.:1.1832
## Max. :285.0 Max. :236.0 Max. :7.267 Max. :1.4323
## geometry
## LINESTRING :10
## epsg:4326 : 0
## +proj=long...: 0
##
##
##
Se agregaron 2 nuevos campos llamados duration y distance donde podemos ver que:
La duración promedio de los viajes es de 5,63 min, mientras que la distancia promedio es de 1,10 km.
El recorrido más corto es de 0,88 km y el más largo de 1,43 km.
El recorrido de menor duración es en 4,02 min y el de mayor es 7,26 min.
Ahora hagamos un mapa con todos los recorridos:
leaflet(ruteo_top10) %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addPolylines(color = "red",
label = paste("Distancia:", round(ruteo_top10$distance, 2), "|", "Duración:", round(ruteo_top10$duration, 2)))
Y ajustemos un poco la estética como por ejemplo:
paleta <- c(low="gold", high= "deeppink4")
leaflet(ruteo_top10) %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addPolylines(color = ~colorNumeric(paleta, ruteo_top10$distance)(distance),
label = paste("Desde", ruteo_top10$src, "hasta", ruteo_top10$dst, "|", "Distancia:", round(ruteo_top10$distance, 2), "km", "|", "Duración:", round(ruteo_top10$duration, 2), "min")) %>%
addLegend("bottomright", pal = colorNumeric(paleta, ruteo_top10$distance), values = ~distance,
title = "Distancia",
labFormat = labelFormat(suffix = "km")) %>%
addCircleMarkers(data=top_10, ~start_x, ~start_y,
popup = paste(ruteo_top10$src, top_10$start_station_name),
color="white",
fillOpacity = 1,
radius=6)%>%
addCircleMarkers(data=top_10, ~end_x, ~end_y,
popup = paste(ruteo_top10$dst, top_10$end_station_name),
color="white",
fillOpacity = 1,
radius=6)
Por último, probemos cambiando el mapa de fondo por uno más claro:
leaflet(ruteo_top10) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolylines(color = ~colorNumeric(paleta, ruteo_top10$distance)(distance),
label = paste("Desde", ruteo_top10$src, "hasta", ruteo_top10$dst, "|", "Distancia:", round(ruteo_top10$distance, 2), "km", "|", "Duración:", round(ruteo_top10$duration, 2), "min")) %>%
addLegend("bottomright", pal = colorNumeric(paleta, ruteo_top10$distance), values = ~distance,
title = "Distancia",
labFormat = labelFormat(suffix = "km")) %>%
addCircleMarkers(data=top_10, ~start_x, ~start_y,
popup = paste(ruteo_top10$src, top_10$start_station_name),
color="gray",
fillOpacity = 1,
radius=6)%>%
addCircleMarkers(data=top_10, ~end_x, ~end_y,
popup = paste(ruteo_top10$dst, top_10$end_station_name),
color="gray",
fillOpacity = 1,
radius=6)
Elegir y descargar de algun portal de datos abiertos un dataset que contenga servicios esenciales (hospitales, escuelas, comisarías, etc) de alguna Ciudad que sea de su interés.
Elegir un barrio cualquiera de la Ciudad en cuestión y calcular el centroide con st_centroid(). Los datos pueden ser descargados del portal open data de la ciudad y deben tener ubicación geográfica.
Estimar la distancia entre el centroide calculado en el punto 2 y los ítems que componen la capa descargada (punto 1). Describir los resultados obtenidos y hacer un mapa con los ruteos. Para resolver este ejercicio pueden ver el siguiente ejemplo.