Distancias medidas en red de calles
Entre grilla de puntos lat-long y tramos de calle, not-chajá mode

* Un chajá o Spiderman, ambos ejemplos sirven, pero prefiero no mencionar marcas.

El objetivo del presente trabajo será obtener las distancias entre cada elemento de una grilla de puntos lat-long y sus calles aledañas, ambos correspondientes a la Ciudad de Buenos Aires. Se usará fuertemente el servicio de ruteo de OSRM, cargado en un servidor local según lo aprendido en este tutorial.

Paso 1: Carga de paquetes y fuentes de datos

Cargamos los paquetes a utilizar: dplyr, sf, here y ggplot2. Se utilizarán otros en momentos específicos y será indicado como corresponde.

library(dplyr)
library(sf)
library(here)
library(ggplot2)

Luego es necesario cargar los dataset seleccionados. La grilla la generaremos desde cero, y sí cargaremos el dataset de calles.

Cargaremos también el dataset de barrios para filtrar puntos de la grilla,

El callejero se puede obtener en BA Data, pero como no quiero descargarlo cada vez que ejecuto lo descargué.

callejero_path <- here("data", "callejero.csv")
calles <- st_read(callejero_path,
                  stringsAsFactors = FALSE)
## Reading layer `callejero' from data source `/home/rstudio/work/data/callejero.csv' using driver `CSV'
## Simple feature collection with 30472 features and 30 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -58.53244 ymin: -34.70574 xmax: -58.34191 ymax: -34.52947
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=-34.6297166 +lon_0=-58.4627 +k=0.9999980000000001 +x_0=100000 +y_0=100000 +ellps=intl +units=m +no_defs
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
# Para etapa de testeo me quedo sólo con Belgrano
barrios <- barrios %>% 
  filter(barrio == "BELGRANO")

rm(callejero_path)

Importar columnas como factor tiene el beneficio de requerir menos espacio, pero puede inducir errores como importar un ID numérico como factor y que lo que se vea no es la operación real. Siempre que mi memoria pueda manejarlo prefiero los string.

Las calles no poseen SRID, los barrios son polígonos con SRID 4326.

Paso 2: Armado de la grilla de puntos x y

Con los barrios de la ciudad obtendré el bounding box y también usaré la geometría de CABA para filtrar espacialmente los puntos de la grilla.

Genero la geometría de la ciudad a partir de los barrios.

caba <- barrios %>%
    mutate(ciudad = "CABA") %>%
    group_by(ciudad) %>%
    summarise() %>%
    smoothr::fill_holes(threshold = 0.01)

ggplot() + 
  geom_sf(data = caba, aes())

Listo lo que queríamos. La función fill_holes() fue usada para eliminar algunos huecos entre barrios que no corresponden, donde falló el snapping de las geometrías de cada barrio.

Obtengo el bounding box de CABA.

bbox <- st_bbox(caba)
bbox %>% round(3)
##    xmin    ymin    xmax    ymax 
## -58.473 -34.575 -58.425 -34.532

Luego genero las secuencias de puntos que conformarán la grilla. La resolución elegida es de latitudes con 3 decimales, lo que implica unos 110 m de separación entre puntos contiguos.

long <- seq(from = bbox[1] %>% round(3),
            to = bbox[3] %>% round(3),
            by = 0.001)

lat <- seq(from = bbox[2] %>% round(3),
           to = bbox[4] %>% round(3),
           by = 0.001)

puntos_grilla <- data.table::CJ(long = long,
                                lat = lat) %>%
  st_as_sf(coords = c("long", "lat"),
           crs = 4326)

ggplot() + 
    geom_sf(data = caba, aes()) +
    geom_sf(data = puntos_grilla,
          aes(),
          size = 0.3,
          alpha = 0.3)

rm(long, lat, bbox)

# library(data.table)
# CJ() "cross join"
# detach("package:data.table")

No cargo el paquete solo para esa función, pues se pisa con algunas funciones de dplyr. Aún no soy muy ducho con el detach() para evitar enmascaramientos posteriores.

Dado que no quiero quedarme con puntos de la bbox que estén por fuera de CABA, los filtro espacialmente.

puntos_grilla <- st_join(puntos_grilla, caba) %>%
    filter(!is.na(ciudad)) %>%
    select(-ciudad) %>%
    mutate(id_grilla = row_number())

# opción de sample para testeo:
# set.seed(94062)
# origenes <- puntos_grilla %>%
#     sample_n(size = 100)

origenes_4326 <- puntos_grilla

origenes_5347 <- puntos_grilla %>% 
  st_transform(5347)

origenes_4326_list <- cbind(origenes_4326,
                            st_coordinates(origenes_4326)) %>%
  st_drop_geometry()

rm(puntos_grilla)

Aquí puede realizarse una muestra de los puntos para testear tiempos antes de tomar todos los datos. En total la grilla de CABA posee 20015 puntos.

También genero un ID de grilla que será el identificador de cada punto.

Vuelvo a chequear grilla filtrada.

ggplot() + 
  geom_sf(data = caba, aes()) +
  geom_sf(data = origenes_4326,
          aes(),
          size = 0.2,
          alpha = 0.3)

Paso 3: Preparación de dataset de calles

Reviso el callejero que cargué.

ggplot() + 
  geom_sf(data = caba,
          aes(),
          size = 0.2,
          alpha = 0.3) +
  geom_sf(data = calles,
          aes())

Hay algo raro. Pruebo nuevamente seteando el CRS.

st_crs(calles) <- 4326
# No comprendo por qué no trae EPSG code, así que lo seteo.
# No es proyección GKBA.

ggplot() + 
  geom_sf(data = caba,
          aes(),
          size = 0.2,
          alpha = 0.3) +
  geom_sf(data = calles,
          aes())

Ahora la geom se ve OK.

Revisemos los atributos. Para esto puede hacerse un head(10) o un str().

calles %>% str()
## Classes 'sf' and 'data.frame':   30472 obs. of  31 variables:
##  $ WKT       : chr  "LINESTRING (-58.4693950261311 -34.5908231628172,-58.4694624545086 -34.590981184668,-58.4696694802765 -34.591727"| __truncated__ "LINESTRING (-58.4292181631792 -34.5508470510858,-58.4296918723879 -34.5505234549563)" "LINESTRING (-58.4111723965244 -34.5576823246023,-58.4121097005867 -34.5570327992133,-58.4127946167632 -34.55660"| __truncated__ "LINESTRING (-58.4843475765456 -34.5576795173653,-58.484638641402 -34.5578139064837)" ...
##  $ id        : chr  "28686" "896" "1724" "1729" ...
##  $ codigo    : chr  "17138" "16003" "16003" "19028" ...
##  $ nomoficial: chr  "PUNTA ARENAS" "OBLIGADO RAFAEL, Av.Costanera" "OBLIGADO RAFAEL, Av.Costanera" "LARRALDE, CRISOLOGO AV." ...
##  $ alt_izqini: chr  "0" "6182" "4202" "4202" ...
##  $ alt_izqfin: chr  "0" "6200" "4500" "4220" ...
##  $ alt_derini: chr  "0" "6181" "4201" "4201" ...
##  $ alt_derfin: chr  "0" "6199" "4499" "4219" ...
##  $ nomanter  : chr  "" "COSTANERA NORTE, Av." "COSTANERA NORTE, Av." "REPUBLIQUETAS" ...
##  $ nom_mapa  : chr  "TÚNEL PUNTA ARENAS" "AV.COSTANERA RAFAEL OBLIGADO" "AV.COSTANERA RAFAEL OBLIGADO" "AV. CRISOLOGO LARRALDE" ...
##  $ tipo_c    : chr  "TÚNEL" "AVENIDA" "AVENIDA" "AVENIDA" ...
##  $ long      : chr  "334.60341680080836" "56.38" "297.26" "30.59" ...
##  $ sentido   : chr  "DOBLE" "DOBLE" "DOBLE" "CRECIENTE" ...
##  $ cod_sent  : chr  "2" "2" "2" "1" ...
##  $ observa   : chr  "Viaducto - Túnel inaugurado en Abril de 2009" "" "" "" ...
##  $ bicisenda : chr  "-" "Ciclovías" "Ciclovías" "-" ...
##  $ lado_ciclo: chr  "" "" "" "" ...
##  $ recorrid_x: chr  "" "" "" "" ...
##  $ ciclo_obse: chr  "" "Construcción Año 2014" "Construcción Año 2014" "" ...
##  $ tooltip_bi: chr  "" "Ciclovía" "Ciclovía" "" ...
##  $ red_jerarq: chr  "VÍA DISTRIBUIDORA COMPLEMENTARIA" "VÍA DISTRIBUIDORA PRINCIPAL" "VÍA DISTRIBUIDORA PRINCIPAL" "VÍA DISTRIBUIDORA PRINCIPAL" ...
##  $ red_tp    : chr  "" "" "" "" ...
##  $ ffcc      : chr  "SI" "" "" "SI" ...
##  $ tipo_ffcc : chr  "Túnel" "" "" "Paso a Nivel" ...
##  $ comuna    : chr  "15" "13" "14" "12" ...
##  $ com_par   : chr  "15" "13" "14" "12" ...
##  $ com_impar : chr  "15" "13" "14" "12" ...
##  $ barrio    : chr  "PATERNAL" "PALERMO" "PALERMO" "SAAVEDRA" ...
##  $ barrio_par: chr  "PATERNAL" "PALERMO" "PALERMO" "SAAVEDRA" ...
##  $ barrio_imp: chr  "PATERNAL" "PALERMO" "PALERMO" "SAAVEDRA" ...
##  $ geometry  :sfc_LINESTRING of length 30472; first list element:  'XY' num [1:6, 1:2] -58.5 -58.5 -58.5 -58.5 -58.5 ...
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "WKT" "id" "codigo" "nomoficial" ...

Me quedo solo con los atributos que utilizaré. Por otra parte, para las operaciones espaciales utilizaré coordenadas planas (CRS 5347), por lo que armo el dataset de calles y centroides en esas coordenadas.

Los centroides son calculados mejor si se aplican sobre coordenadas planas, también.

calles <- calles %>%
    select(id, geometry, comuna, barrio) %>%
    mutate(id = id %>% as.integer()) %>%
    rename(id_calle = id)

destinos_5347 <- calles %>%
  st_transform(crs = 5347) %>% 
  st_centroid()
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
destinos_4326 <- destinos_5347 %>%
  st_transform(crs = 4326)
# con esto busco hacer centroides siempre sobre coordenadas planas, no con datos en EPSG 4326.

destinos_4326_list <- cbind(destinos_4326,
                       st_coordinates(destinos_4326)) %>%
    st_drop_geometry()

rm(calles)

Hubo un problema con el orden de operaciones, y pasar a 5347 y volver a 4326 me tiraba los puntos en cualquier lado.

Veo efectivamente que los puntos están dentro de CABA.

st_crs(origenes_4326)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(origenes_5347)
## Coordinate Reference System:
##   EPSG: 5347 
##   proj4string: "+proj=tmerc +lat_0=-90 +lon_0=-60 +k=1 +x_0=5500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
st_crs(destinos_4326)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(destinos_5347)
## Coordinate Reference System:
##   EPSG: 5347 
##   proj4string: "+proj=tmerc +lat_0=-90 +lon_0=-60 +k=1 +x_0=5500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
ggplot() + 
  geom_sf(data = origenes_5347 %>%
            st_buffer(dist = 600) %>% 
            head(1000),
          aes(),
          alpha = 0.5) +
  geom_sf(data = destinos_5347 %>%
            head(1000),
          aes(),
          alpha = 0.5)

Chequeo más feo ever, pero makes sense y sigo hacia adelante.

Con esto tengo armado mis datos orígenes y destinos. Puntos importantes: para ruteo OSRM requiero pasar latitud y longitud, (EPSG 4326), mientras que para operaciones geoespaciales es deseable utilizar coordenadas planas (EPSG 5347, POSGAR 2007).

Paso 4: Armado de puntos de destinos para cada elemento de punto grilla

Para cada punto de la grilla, realizaré el ruteo a todos los callejeros cuyo centroide esté en un radio de 600 metros. Para ello, primero hago buffers de puntos grilla y los interseco con los centroides de calles.

t_ini <- Sys.time()

interseccion_od <- st_intersection(origenes_5347 %>%
                                     st_buffer(dist = 600) %>% 
                                     select(id_grilla, geometry),
                                   destinos_5347 %>% 
                                     select(id_calle, geometry)) %>% 
    st_drop_geometry()
# intesección buffer con punto: punto.
# hago drop geometry pues no usaré el punto que se generó.

interseccion_od %>% nrow()
## [1] 106378
t_fin <- Sys.time()

print(t_fin-t_ini)
## Time difference of 8.128757 secs

Al seleccionar solo id de cada tabla (y la geometría que lo hace por defecto), la tabla de intersección queda más comprensible.

Es esta tabla la que tendra una fila por ruteo que haga, es decir que serán en total 3.272.212 llamados al servidor (toda la ciudad, para Belgrano serán 106.378).

Paso 5: Ruteo entre par de puntos

Cargo el paquete osrm y lo configuro para pegarle a mi servidor local.

También genero el formato de tabla output que preciso.

Luché un rato con la IP local, una opción es armar la URL con localhost, o desde Ubuntu colocar en la terminal hostname -I para ver la IP que tienen. También desde la terminal, ifconfig me mostró esa IP bajo el un adaptador llamado wlp2s0. Estando WiFi variables entiendo que corresponde chequear cada vez o ir por la URL genérica del localhost.

En el siguiente chuck estuve un rato probando URLs y con un caso de juguete, mientras veía como reaccionaba la terminal donde había levantado el servidor. Estuve un montón hasta darme cuenta que la URL del parámetro osrm.server tiene que venir con una / al final.

library(osrm)

# options(osrm.server = "http://10.78.17.102:5000/")
# options(osrm.server = "http://127.0.0.1:5000/")
# options(osrm.server = "http://localhost:5000/")
# options(osrm.server = "http://192.168.0.107:5000/")
options(osrm.server = "http://172.17.0.1:5000/")

ciudades <- read.csv("http://bitsandbricks.github.io/data/aglomerados.csv")

distancias <- osrmTable(loc = ciudades)

ruta <- osrmRoute(src = c(punto = "hogar_dulce_obelisco",
                          lon = -58.3824556,
                          lat = -34.6050104),
                  dst = c(punto = "el_rio",
                          lon = -58.4675927,
                          lat = -34.5247167), 
                  returnclass = "sf",
                  overview = "full")

ggplot() +
  geom_sf(data = st_read("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios.geojson")) +
  geom_sf(data = ruta,
          size = 1.5,
          alpha = 0.7)
## 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

rm(ruta, distancias, ciudades)

En búsqueda de la optimización, decidí unir de antemano datos espaciales con los que llamaré al servidor.

# uno campos XY origen
interseccion_od <- left_join(interseccion_od,
                             origenes_4326_list,
                             by = c("id_grilla" = "id_grilla"))

# uno campos XY destino
interseccion_od <- left_join(interseccion_od,
                             destinos_4326_list %>%
                                 select(id_calle, X, Y),
                             by = c("id_calle" = "id_calle"),
                             suffix = c("_grilla", "_calle"))

# lo uno ahora pues intersección es en 5347 y yo busco lat-long

En este momento venía el momento poco ortodoxo de un loop en R. Sin embargo la autoregulación de las circunstancias y su mala performance me llevaron a replicar el método de este ejemplo, compartilhado pela minha companheira Angie.

ruteo_grilla_calle_fun <- function(id_grilla_o,
                                   o_x, o_y,
                                   id_calle_d,
                                   d_x, d_y) {
  
  ruta <- osrmRoute(src = c(id_grilla_o, o_x, o_y),
                    dst = c(id_calle_d, d_x, d_y), 
                    returnclass = "sf",
                    overview = "full")
  
  return(cbind(id_grilla_origen = id_grilla_o,
               id_calle_destino = id_calle_d,
               ruta))
}

No comprendo por qué los parámetros src y dst requieren un nombre que luego vuelvo a cbindear.

Con el siguiente código realizo el ruteo.

Este es mi cuello de botella que aún no tengo resuelto. ¿pmap? ¿reduce?

ruteo_grilla_calle <- list(interseccion_od$id_grilla,
                           interseccion_od$X_grilla,
                           interseccion_od$Y_grilla,
                           interseccion_od$id_calle,
                           interseccion_od$X_calle,
                           interseccion_od$Y_calle)

# testing case:
# me quedo con 1000 filas
ruteo_grilla_calle <- lapply(X = ruteo_grilla_calle,
                             FUN = head,
                             n = 1000)

t_ini <- Sys.time()

output <- purrr::pmap(.l = ruteo_grilla_calle,
                      .f = ruteo_grilla_calle_fun)

t_fin <- Sys.time()

print("tiempo de ruteo:")
## [1] "tiempo de ruteo:"
print(t_fin-t_ini)
## Time difference of 8.592176 secs
t_ini <- Sys.time()

output <- purrr::reduce(.x = output,
                        .f = bind_rows)
# output es de tipo sf. that's so nice.

# pmap_dfr debería recibir los mismos inputs y ahorrarme el reduce, pero no sé si preserva el tipo SF

# output <- purrr::pmap_dfr(ruteo_grilla_calle,
#                           ruteo_grilla_calle_fun)

t_fin <- Sys.time()

print("tiempo de reduce:")
## [1] "tiempo de reduce:"
print(t_fin-t_ini)
## Time difference of 0.2965205 secs

La función pmap toma 2 argumentos: .l es una lista de vectores, y .f la función a aplicar, tomando como parámetros los n-ésimos elementos de cada vector de la lista.

Tomando un head(15000) y reduciendo con rbind lo interrumpí a los 38 minutos por extremandamente lento. Lo mejoré un poco al realizar la reducción con bind_rows. Go team dplyr::bind_rows vs. team base::rbind.

Mi conocimiento limitado de R me dice que pmap paraleliza en los cores de mi PC, por lo que pude ver haciendo %sh htop. Nuevamente viendo htop, intuyo que el reduce corre en un único core que se pone al 100%.

Tomando un head(1000) corre en unos 9.3 s, not bad para testear. Todo Belgrano (106.000 ruteos) da unos tiempo_estimado_eterno minutos. Nota: cuando crece n el reduce demora mucho más que el pmap y lo que pensé que serían 20 minutos fue algo que interrumpí por eterno. Not nice.

Pasar de 106.000 a los 3.272.212 daría un tiempo total estimado de tiempo_eterno_en_horas horas.

Sólo para la posteridad y el hall of shame de past-Andrew, comparto el loop que intenté y demoraba más de lo deseado. Releer old shitty code sirve para subir la confianza y conscientizar las competencias.

t_ini <- Sys.time()

duracion_v <- numeric()
distancia_red_v <- numeric()
ruta_od_v <- character()

n_total <- nrow(interseccion_od)

for (od in 1:n_total) {
  
  id_grilla <- interseccion_od[od, ]$id_grilla
  locate_bool_o <- origenes_4326_list$id_grilla == id_grilla
  p_o_x <- origenes_4326_list[locate_bool_o, "X"]
  p_o_y <- origenes_4326_list[locate_bool_o, "Y"]

  id_calle <- interseccion_od[od, ]$id_calle
  locate_bool_d <- destinos_4326_list$id_calle == id_calle
  p_d_x <- destinos_4326_list[locate_bool_d, "X"]
  p_d_y <- destinos_4326_list[locate_bool_d, "Y"]

  ruta_od <- osrmRoute(src = c(p_o_x, p_o_y),
                       dst = c(p_d_x, p_d_y),
                       returnclass = "sf",
                       overview = "full")
  
  duracion_v <- c(duracion_v, ruta_od$duration)
  distancia_red_v <- c(distancia_red_v, ruta_od$distance)
  ruta_od_v <- c(ruta_od_v, ruta_od$geometry)
  
  print(od)
  print(paste0(od/n_total*100, "%"))
  
}
# osrm devuelve info distinta que el ruteador de USIG

output <- data.frame(id_grilla = interseccion_od$id_grilla,
                     id_calle = interseccion_od$id_calle,
                     duracion = duracion_v,
                     distancia_red = distancia_red_v,
                     ruta_od = ruta_od_v %>% as.character())

t_fin <- Sys.time()

print(t_fin - t_ini) 

rm(id_grilla, id_calle, locate_bool_o, locate_bool_d, p_o_x, p_o_y, p_d_x, p_d_y, od, duracion_v, distancia_red_v, ruta_od_v, t_ini, t_fin)

Cuando lo corrí, lo interrumpí antes de que terminara, pues hizo corrió cerca de 50.000 requests en más de 11 minutos. Estimé un proceso total de 12 horas en correr los 3.272.212 de llamados, y me dije “tengo que optimizarlo”. La optimización implicó su destrucción y reemplazo por el pmap.

El estimado HAVB en su post habla de 15K ruteos por minuto, es decir que yo debería bajar a menos de 4 horas, o al menos tener un orden similar. Trataré de compensar su expertise con mi hardware (?).

3272212/15000/60
## [1] 3.635791
# total de filas / filas por minuto / minutos por hora

Nota de Andy: Estaría bueno ver una barra incrementándose para bajar la incertidumbre de tiempo.

Finalmente, viene el testeo, donde agrego las geometrías de OD + ruta.

Aunque tenga la geometría de ruta, las rutas se superponen, por lo que las geoms de OD son útiles.

t_ini <- Sys.time()

# pensaba añadir las geometrías pero voy por columnas lat-long

output <- left_join(x = interseccion_od,
                    y = output,
                    by = c("id_grilla" = "id_grilla_origen",
                           "id_calle" = "id_calle_destino"))
# esto debería unirme duration, distance, geometry
# no comprendo por qué no me deja hacer: 
# output %>% select(-c("src", "dst"))
# Error: Only strings can be converted to symbols. --> fuck this

output <- output %>%
    rename(distancia_red = distance)

# uno geom origen (en coordenadas planas)
# output <- left_join(output,
#                     origenes_5347 %>% as.data.frame(),
#                     by = c("id_grilla_origen" = "id_grilla"),
#                     suffix = c("_ruta", "_grilla")
#                     )

# uno geom destino (en coordenadas planas)
# output <- left_join(output,
#                     destinos_5347 %>% 
#                       select(id_calle, geometry, barrio, comuna),
#                     by = c("id_calle" = "id_calle")
#                     )

# output <- output %>% 
#     rename(geometry_calle = geometry)
# genero distancia a vuelo de pájaro
output <- output %>%
    mutate(distancia_euclideana = st_distance(lapply(X = output %>% 
                                                         pull(id_grilla),
                                                     FUN = function(x) origenes_5347[id_grilla == x, "geometry"]),
                                              lapply(X = output %>% 
                                                         pull(id_calle),
                                                     FUN = function(x) destinos_5347[id_calle == x, "geometry"]),
                                              by_element = TRUE) %>% as.numeric())
# validar que el lapply funciona OK, no sé qué onda con las geometries.
# eterno para 1000 valores, revisar.

output <- output %>%
    mutate(geom_grilla = st_point(x = c(X_grilla, Y_grilla),
                                  dim = "XY",
                                  crs = 4326),
           geom_calle = st_point(x = c(X_calle, Y_calle),
                                  dim = "XY",
                                  crs = 4326))
# llevo distancia a metros
output <- output %>%
  mutate(distancia_red = distancia_red*1000)

Hago primera visualización.

WIP, aun está parte está en producción porque no superé mi cuello de botella de arriba.

library(leaflet)

data <- output %>% 
  filter(distancia_red <= 1000) %>% 
  filter(id_grilla_origen %in%
           sample(x = output %>% 
                    pull(id_grilla_origen),
                  size = 15))

leaflet(data = data,
        width = "100%") %>%
  addTiles() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolylines(color = ~colorNumeric("viridis",
                                     data$distancia_red)(distancia_red),
               popup = data$distancia_red %>%
                   lapply(htmltools::HTML),
               label = ~(paste0(id_grilla, "->", id_calle) %>% as.character())
               opacity = 0.3,
               group = "Rutas") %>%
  addLegend("bottomright",
            pal = colorNumeric("viridis",
                               data$distancia_red),
            values = ~distancia_red,
            title = "Distancia",
            opacity = 0.8) %>% 
  addCircleMarkers(lng = ~X_grilla,
                   lat = ~Y_grilla,
                   radius = 1/distancia_red,
                   color = ~colorNumeric("viridis",
                                     data$distancia_red)(distancia_red),
                   fillOpacity = 0.8,
                   label = ~(id_grilla %>% as.character()),
                   group = "Grilla origen") %>%
  addLegend("bottomright", pal = pal2, values = ~en_cabecera, title = "En cabecera", opacity = 1) %>% 


# necesito agregar puntos (x, y) de grilla origen

Visualizo la relación entre la distancia de ruteo y la distancia por aire.

output <- output %>% 
    mutate(dist_ratio = distancia_euclideana/distancia_red)

p <- ggplot(data = output) +
    geom_point(aes(x = distancia_red,
                   y = distancia_euclideana),
               size = 0.8,
               alpha = 0.3)
p

q <- ggplot(data = output) +
    geom_point(aes(x = distancia_red,
                   y = distancia_euclideana),
               size = 0.8,
               alpha = 0.3) +
    facet_wrap(. ~ barrio)
q

r <- ggplot(data = output) +
    geom_boxplot(aes(y = dist_ratio,
                     group = barrio))
r

s <- ggplot(data = output) +
    geom_histogram(aes(x = dist_ratio),
                   bins = 30) +
    facet_wrap(. ~ barrio)
s

rm(p, q, r, s)

Finalmente, escribo el output en un archivo CSV y en un objeto .RData.

# revisar en este punto qué columnas quiero exportar.
# tal vez exportar diccionarios OD en lugar de columnas unidaes.

write.csv2(x = output,
           file = here("data", "output_dist_ruteo.csv"),
           row.names = FALSE,
           fileEncoding = "UTF-8",
           sep = ";")

save(output,
     file = here("data", "output_dist_ruteo.RData"))
Observaciones:
Pendiente