library("tidyverse")
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library("ggplot2")
library("sf")
## Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
library("osmdata")
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library("rmarkdown")
library("lubridate")
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library("ggmap") 
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library("mapview")
library("hrbrthemes")
## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
##       Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
##       if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow
library("RColorBrewer")
library("leaflet")
library("osmextract")
## Warning: package 'osmextract' was built under R version 4.2.1
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library("forcats")
library("osrm")
## Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
## Routing: OSRM - http://project-osrm.org/
## sp support will be dropped in the next major release, please use sf objects instead.
library("hereR")
## Warning: package 'hereR' was built under R version 4.2.2
## 
## Attaching package: 'hereR'
## The following objects are masked from 'package:ggmap':
## 
##     geocode, route
Nordelta = data.frame(-58.6814666,-34.4000434)

colnames(Nordelta) <- c('lon','lat')

Nordelta =  st_as_sf(x = Nordelta, coords = c("lon", "lat"),  crs = 4326)
Recoleta =  data.frame(-58.3914178,-34.5911356)

colnames(Recoleta) <- c('lon','lat')

Recoleta =  st_as_sf(x = Recoleta, coords = c("lon", "lat"),  crs = 4326)
mapview(Recoleta) + mapview(Nordelta)
Nordelta_Isocrona = osrmIsochrone(loc = Nordelta, seq(from = 0, to = 20, seq = 5),  osrm.profile = "foot")
## Warning: In seq.default(from = 0, to = 20, seq = 5) :
##  extra argument 'seq' will be disregarded
Recoleta_Isocrona = osrmIsochrone(loc = Recoleta, seq(from = 0, to = 20, seq = 5),  osrm.profile = "foot")
## Warning: In seq.default(from = 0, to = 20, seq = 5) :
##  extra argument 'seq' will be disregarded
mapview(Nordelta, color = "red") + mapview(Recoleta_Isocrona, zcol = "center") + mapview(Nordelta_Isocrona,zcol = "center") + mapview(Recoleta, color = "red")
st_area(Nordelta_Isocrona)
## Units: [m^2]
##  [1]   1057.182   7418.705  12441.657  13567.240  18559.638  24614.550
##  [7]  15186.819  41025.561  31143.952  58031.666  50053.955  65034.147
## [13]  85974.506 164094.656 164199.619 208573.088 229244.279 225769.725
## [19] 275322.458 343645.741
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
st_area(Recoleta_Isocrona)
## Units: [m^2]
##  [1]   1377.216  56952.430  64529.277  88673.466 154156.213 151480.890
##  [7] 155590.650 157078.170 207506.255 207939.890 240066.494 267183.166
## [13] 287232.381 279834.471 280230.254 329130.401 333451.935 334838.506
## [19] 405701.162

Nos quedamos con el más grande de cada dataset.

Nordelta_Isocrona_top = osrmIsochrone(loc = Nordelta, seq(from = 0, to = 20),  osrm.profile = "foot")
Recoleta_Isocrona_top = osrmIsochrone(loc = Recoleta, seq(from = 0, to = 20),  osrm.profile = "foot")
mapview(Nordelta, color = "red") + mapview(Recoleta_Isocrona_top) + mapview(Nordelta_Isocrona_top) + mapview(Recoleta, color = "red")
mapview(Recoleta_Isocrona_top, zcol = "center")
mapview(Nordelta_Isocrona_top, zcol = "center")
st_area(Nordelta_Isocrona_top)
## Units: [m^2]
##  [1]   1056.967   7417.200  12439.133  13564.487  18555.871  24609.552
##  [7]  15183.736  41017.243  31137.633  58019.877  50043.793  65020.942
## [13]  85957.049 164061.328 164166.325 208530.768 229197.814 225723.954
## [19] 275266.647 343576.019
st_area(Recoleta_Isocrona_top)
## Units: [m^2]
##  [1]   1377.216  56952.430  64529.277  88673.466 154156.213 151480.890
##  [7] 155590.650 157078.170 207506.255 207939.890 240066.494 267183.166
## [13] 287232.381 279834.471 280230.254 329130.401 333451.935 334838.506
## [19] 405701.162
st_area(Nordelta_Isocrona_top) / st_area(Recoleta_Isocrona_top)
## Warning in `/.default`(st_area(Nordelta_Isocrona_top),
## st_area(Recoleta_Isocrona_top)): longitud de objeto mayor no es múltiplo de la
## longitud de uno menor
## Units: [1]
##  [1]   0.76746681   0.13023501   0.19276727   0.15297121   0.12037057
##  [6]   0.16245978   0.09758771   0.26112631   0.15005636   0.27902235
## [11]   0.20845805   0.24335718   0.29925960   0.58627991   0.58582656
## [16]   0.63358100   0.68734888   0.67412783   0.67849608 249.47147031
Nordelta_Isocrona_top_auto = osrmIsochrone(loc = Nordelta, seq(from = 0, to = 15),  osrm.profile = "car")
Recoleta_Isocrona_top_auto = osrmIsochrone(loc = Recoleta, seq(from = 0, to = 15),  osrm.profile = "car")
mapview(Nordelta, color = "red") + mapview(Nordelta_Isocrona_top_auto, zcol = "center") 
mapview(Recoleta_Isocrona_top_auto,zcol = "center") + mapview(Recoleta, color = "red")
st_area(Nordelta_Isocrona_top_auto) / st_area(Recoleta_Isocrona_top_auto)
## Warning in `/.default`(st_area(Nordelta_Isocrona_top_auto),
## st_area(Recoleta_Isocrona_top_auto)): longitud de objeto mayor no es múltiplo de
## la longitud de uno menor
## Units: [1]
##  [1] 0.001753565 0.157896948 0.559841197 0.304828806 0.456478488 0.340373182
##  [7] 0.667038384 0.551624338 0.875924544 0.856370303 0.659361358 7.868394917
## [13] 6.810453968

Ahora que ya jugamos un poco con los mapas y las isocronas vamos a ver como varía el acceso a diferentes bienes y servicios según el barrio de CABA.

Vamos a trabajar con la infraestructura universitaria y los barrios

barrios <- st_read('https://bitsandbricks.github.io/data/CABA_barrios.geojson')
## Reading layer `CABA_barrios' from data source 
##   `https://bitsandbricks.github.io/data/CABA_barrios.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 48 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53152 ymin: -34.70529 xmax: -58.33514 ymax: -34.52754
## Geodetic CRS:  WGS 84

Cambiamos el nombre de villa lugano

barrios$BARRIO[barrios$BARRIO == "VILLA LUGANO"] <- 'Villa Lugano'
barrios$BARRIO[barrios$BARRIO == "RECOLETA"] <- 'Recoleta'

Calculamos los centroides de los barrios

barrios_centroides = st_centroid(barrios)
## Warning in st_centroid.sf(barrios): st_centroid assumes attributes are constant
## over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for longitude/
## latitude data

Descomprimimos la geometría

barrios_centroides <- barrios_centroides %>%
mutate(lat = unlist(map(barrios_centroides$geometry,1)),
long = unlist(map(barrios_centroides$geometry,2)))
Lugano = filter(barrios_centroides, BARRIO == "Villa Lugano")
Recoleta = filter(barrios_centroides, BARRIO == "Recoleta")

Universidades

Universidades = read.csv("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/ministerio-de-educacion/universidades/universidades.csv")
colnames(Universidades)[16] <- "BARRIO"
uba = filter(Universidades, universida =="Universidad de Buenos Aires")
library("hereR")
uba = uba %>% st_as_sf(coords = c("long", "lat"), crs = 4326)

Nos pide que sean del mismo tamaño, multiplicamos el dataset de Recoleta

Recoleta = rbind(Recoleta, Recoleta[rep(1, 24), ])
route(Recoleta, uba)
## Simple feature collection with 25 features and 12 fields
## Geometry type: LINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: -58.5235 ymin: -34.6803 xmax: -58.36849 ymax: -34.53432
## z_range:       zmin: 21 zmax: 50
## Geodetic CRS:  WGS 84
## First 10 features:
##    id rank section           departure             arrival    type mode
## 1   1    1       1 2022-11-29 09:34:51 2022-11-29 09:47:25 vehicle  car
## 2   2    1       1 2022-11-29 09:34:51 2022-11-29 09:40:44 vehicle  car
## 3   3    1       1 2022-11-29 09:34:51 2022-11-29 09:55:54 vehicle  car
## 4   4    1       1 2022-11-29 09:34:51 2022-11-29 09:40:15 vehicle  car
## 5   5    1       1 2022-11-29 09:34:51 2022-11-29 09:37:59 vehicle  car
## 6   6    1       1 2022-11-29 09:34:51 2022-11-29 09:40:39 vehicle  car
## 7   7    1       1 2022-11-29 09:34:51 2022-11-29 09:57:30 vehicle  car
## 8   8    1       1 2022-11-29 09:34:51 2022-11-29 09:49:38 vehicle  car
## 9   9    1       1 2022-11-29 09:34:51 2022-11-29 09:47:22 vehicle  car
## 10 10    1       1 2022-11-29 09:34:51 2022-11-29 09:48:49 vehicle  car
##    distance duration duration_base consumption tolls
## 1      3461      754           543      2.1449     0
## 2      2027      353           313      1.1399     0
## 3      5586     1263           744      3.6647     0
## 4      1897      324           284      1.0649     0
## 5       866      188           132      0.4733     0
## 6      1954      348           306      1.1225     0
## 7      8588     1359           922      4.8839     0
## 8      7770      887           706      3.5264     0
## 9      3873      751           483      2.4017     0
## 10     4463      838           516      2.8309     0
##                          geometry
## 1  LINESTRING Z (-58.39517 -34...
## 2  LINESTRING Z (-58.39517 -34...
## 3  LINESTRING Z (-58.39517 -34...
## 4  LINESTRING Z (-58.39517 -34...
## 5  LINESTRING Z (-58.39517 -34...
## 6  LINESTRING Z (-58.39517 -34...
## 7  LINESTRING Z (-58.39517 -34...
## 8  LINESTRING Z (-58.39517 -34...
## 9  LINESTRING Z (-58.39517 -34...
## 10 LINESTRING Z (-58.39517 -34...
Ruteo = route(Recoleta, uba)

Dividimos la duración que está en segundos para calcular en minutos el tiempo

Ruteo = Ruteo %>% mutate(duration = duration/60)
mapview(Ruteo, zcol = "duration") + mapview(uba)

Vamos con un Ruteo en Bicicleta

Ruteo_tp = route(Recoleta, uba, transport_mode = "bicycle")
## Warning in .parse_response(i, out$responses()[[i]]): https://router.hereapi.com/v8/routes: Request 'id = 6' failed. 
##   Status 429; Too Many Requests; The user has sent too many requests in a given amount of time ("rate limiting") (RFC 6585).
## Warning in .parse_response(i, out$responses()[[i]]): https://router.hereapi.com/v8/routes: Request 'id = 7' failed. 
##   Status 429; Too Many Requests; The user has sent too many requests in a given amount of time ("rate limiting") (RFC 6585).
## Warning in .parse_response(i, out$responses()[[i]]): https://router.hereapi.com/v8/routes: Request 'id = 10' failed. 
##   Status 429; Too Many Requests; The user has sent too many requests in a given amount of time ("rate limiting") (RFC 6585).

Dividimos la duración que está en segundos para calcular en minutos el tiempo

Ruteo_tp = Ruteo_tp %>% mutate(duration = duration/60)
mapview(Ruteo_tp, zcol = "duration") + mapview(uba)

Vamos con Lugano

Lugano = rbind(Lugano, Lugano[rep(1, 24), ])
Ruteo_Lugano = route(Lugano, uba)

Dividimos la duración que está en segundos para calcular en minutos el tiempo

Ruteo_Lugano = Ruteo_Lugano %>% mutate(duration = duration/60)
mapview(Ruteo_Lugano, zcol = "duration") + mapview(uba)
sum(Ruteo_Lugano$duration)
## [1] 719.0667
sum(Ruteo$duration)
## [1] 403.8333
"Mientras que desde el centroide de Villa Lugano se demoran"
## [1] "Mientras que desde el centroide de Villa Lugano se demoran"
print(sum(Ruteo_Lugano$duration)) 
## [1] 719.0667
"minutos, como sumatorio de tiempo de viaje a las sedes de la UBA, desde recoleta estos son"
## [1] "minutos, como sumatorio de tiempo de viaje a las sedes de la UBA, desde recoleta estos son"
print(sum(Ruteo$duration)) 
## [1] 403.8333
"minutos"
## [1] "minutos"
ruteo_universidades <- function(o_nombre, o_x, o_y, d_nombre, d_x, d_y) {
  ruta <- osrmRoute(src = c(o_nombre, o_x, o_y),
                    dst = c(d_nombre, d_x, d_y), 
                    returnclass = "sf",
                    overview = "full")
  
  cbind(ORIGEN = o_nombre, DESTINO = d_nombre, ruta)
}