shortest_route_distance.R

# @jrosell: Calculate shortest distance to fuel station from all interest points (useful for feature engineering!)
# @jrosell: Use setNames to be sure distance works.
# @jrosell: Use rowwise() |> group_split() |> map() |> list_rbind() to be able to get ETA with .progress = TRUE
# Original source by @milos-agathon: https://www.youtube.com/watch?v=jmImC0v6qmU

# ==============================================================================
# Preparations

libs <- c(
  "tidyverse", # Data wrangling and visualization
  "osmdata",   # PDI & roads
  "dodgr",     # Compute shortest route
  "sf",        # Map Shapes
  "maptiles",  # Map  Tiles
  "tidyterra",  # Plot maps
  "glue"
)
installed_libs <- libs %in% rownames(installed.packages())
if (any(installed_libs == FALSE)) {
  install.packages(libs[!installed_libs], dependencies = TRUE)
}
invisible(lapply(libs, library, character.only = TRUE))

# ==============================================================================
# Helper functions

extract_amenity_osm_feature <- \(
    xmin = NULL, ymin = NULL, xmax = NULL, ymax = NULL,
    osm_feature = osmdata::available_tags("amenity")$Value
) {
  stopifnot(all(!is.null(c(xmin, ymin, xmax, max))))
  amenities <- c(
      xmin = xmin, ymin = ymin,
      xmax = xmax, ymax = ymax
    ) |>
    sf::st_bbox(crs = 4326) |>
    sf::st_as_sfc(crs = 4326) |> 
    osmdata::opq(  
      timeout = 180,
      memsize = 104857600
    ) |>
    osmdata::add_osm_feature(
        key = "amenity",
        value = osm_feature # "fuel"
    ) |>
    osmdata::osmdata_sf() 
  do.call(rbind, amenities[c("osm_points")]) |>
    dplyr::select("osm_id", "geometry")
}


calculate_shortest_route_vec <- \(
  df, g , to, .progress = FALSE
) {
  seq_len(nrow(df)) |> 
    purrr::map_dbl(.progress = .progress, \(i) {
        calculate_shortest_route(df[i, ], g, to_fuel)$distance  
    })
}

calculate_shortest_route <- \(
    xy, g , to
) {    
  suppressWarnings({
      route_distances <- dodgr::dodgr_dists(
        graph = g,
        from = xy,
        to = to,
        shortest = T
      )
  })    
  route_distances_df <- route_distances |>
    as.matrix() |>
    as.data.frame() |>
    t()
  shortest_route_distance_df <- cbind(
        edge_id = rownames(
            route_distances_df
        ),
        route_distances_df
    ) |>
    as_tibble() |>
    setNames(
        c("edge_id", "distance")
    ) |>
    dplyr::mutate(
        distance = as.numeric(distance)
    ) |>
    dplyr::arrange(distance) |>
    dplyr::distinct(
        edge_id, distance,
        .keep_all = T
    ) |>
    dplyr::slice_head(n = 1)

  shortest_route_distance_df
}

extract_motorcar_highway_weight <- \(  
    xmin = NULL, ymin = NULL, xmax = NULL, ymax = NULL
) {
  g <- dodgr::weight_streetnet(
    extract_hways(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax),
    wt_profile = "motorcar",
    type_col = "highway"
  )
  g
}


extract_hways <- \(
  xmin = NULL, ymin = NULL, xmax = NULL, ymax = NULL
) {
  stopifnot(all(!is.null(c(xmin, ymin, xmax, max))))
  the_bbox <- c(
      xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax
    ) |>
    sf::st_bbox(crs = 4326) |> # WGS84
    sf::st_as_sfc(crs = 4326)
  the_roads <- osmdata::opq(
        bbox = the_bbox,
        timeout = 180,
        memsize = 104857600
    ) |>
    osmdata::add_osm_feature(
        key = "highway"
    ) |>
    osmdata::osmdata_sf()

  the_roads$osm_lines
}

extract_paths_sf <- \(paths) {
  paths_sf <- lapply(
    paths, function(x) {
        lapply(
            x, function(y) {
                path_xy <- v[match(
                    y, v$id
                ), ]
                sf::st_linestring(
                    as.matrix(
                        path_xy[, c("x", "y")]
                    )
                ) |>
                    sf::st_sfc(crs = 4326)
            }
        )
    }
  )

  paths_sf <- lapply(
    paths_sf, function(x) {
        do.call(
            rbind, x
        )
    }
  )

  paths_sf <- do.call(rbind, paths_sf)

  paths_sf <- sf::st_sf(
    from = rep(
        names(paths),
        each = nrow(xy)
    ),
    to = rep(
        names(paths),
        times = nrow(xy)
    ),
    geometry = paths_sf[, 1],
    crs = 4326
  )
  paths_sf
}

extract_shortest_path <- \(g, xy, to) {
  paths <- dodgr::dodgr_paths(
    graph = g,
    from = xy,
    to = to
  )
  shortest_route_distance_df <- calculate_shortest_route(
    xy, g , to
  )
  paths_sf <- extract_paths_sf(paths)
  df_graph_edges <- paths |>
      unlist() |>
      as.matrix() |>
      as.data.frame()
  df_graph_edges <- cbind(
      id = rownames(
          df_graph_edges
      ), df_graph_edges
  )
  names(df_graph_edges)[2] <- "edge_id"
  shortest_line <- df_graph_edges |>
      dplyr::filter(
          stringr::str_detect(
              id,
              shortest_route_distance_df$edge_id
          )
      )
  shortest_path_sf <- g |> 
    dodgr::dodgr_to_sf() |>
    dplyr::filter(
        to_id %in% unique(
            shortest_line$edge_id
        )
    ) |>
    sf::st_intersection(paths_sf)
  shortest_path_sf
}

extract_tiles <- \(xmin, ymin, xmax, ymax) {
  the_tiles <- c(
      xmin = xmin, ymin = ymin,
      xmax = xmax, ymax = ymax
    ) |>
    sf::st_bbox(crs = 4326) |>
    sf::st_as_sfc(crs = 4326) |> 
    maptiles::get_tiles(
      provider = "CartoDB.Positron",
      zoom = 14,
      crop = T,
      # project = F,
      forceDownload = T
    )
  the_tiles
}


# ==============================================================================
# Example: Plot the shortest path from a random point to the nearest
# fuel station in Belgrade.

xmin <- 20.417487
ymin <- 44.800668
xmax <- 20.484438
ymax <- 44.832196

to_fuel <-
  extract_amenity_osm_feature(
    xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax,
    osm_feature = "fuel"
  ) |>
  sf::st_coordinates()

g <- extract_motorcar_highway_weight(
  xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax
)
v <- dodgr::dodgr_vertices(g)
xy <- v[sample(
  nrow(v),
  size = 1
),]

ggplot() +
  tidyterra::geom_spatraster_rgb(
      data = extract_tiles(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax)
  ) +
  geom_sf(
      data = extract_amenity_osm_feature(
        xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax,
        osm_feature = "fuel"
      ),
      color = "blue",
      inherit.aes = F
  ) +
  geom_sf(
      data = extract_hways(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax),
      color = "black",
      size = .15,
      alpha = .5,
      inherit.aes = F
  ) +    
  geom_sf(
      data = extract_shortest_path(g, xy, to_fuel),
      color = "red",
      size = .5,
      inherit.aes = F
  ) +
  theme_void()

  

# ==============================================================================
# Example: Add to each row the distance to the nearest fuel station.

xmin <- 20.417487
ymin <- 44.800668
xmax <- 20.484438
ymax <- 44.832196

to_fuel <-
  extract_amenity_osm_feature(
    xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax,
    osm_feature = "fuel"
  ) |>
  sf::st_coordinates()
g <- extract_motorcar_highway_weight(
  xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax
)
v <- dodgr::dodgr_vertices(g)




# Attempt 1: 100 in 5.4s, no ETA
tictoc::tic()
results <- v[1:100,] |> 
  nest_by(id) |> 
  mutate(
    distance = calculate_shortest_route(
      data, g, to_fuel      
    )$distance
  ) |>
  unnest(cols = c(data)) |> 
  print()
tictoc::toc()

# Attempt 2: 100 in 5.5s, no ETA
tictoc::tic()
results <- v[1:100,] |> 
  rowwise(id) |> 
  group_map(
    \(df, id) {
      df$distance = calculate_shortest_route(
          df, g, to_fuel
      )$distance
      df
    }    
  ) |> 
  list_rbind() |> 
  print()
tictoc::toc()

# Attempt 3: 100 in 5.5s with ETA
tictoc::tic()
results <- v[1:100,] |> 
  rowwise(id) |> 
  group_split() |> 
  map(
    .progress = TRUE,
    \(df) {
      df$distance = calculate_shortest_route(
            df, g, to_fuel
      )$distance
      df
    }
  ) |> 
  list_rbind() |> 
  print()
tictoc::toc()

# Full calculation: 13447 in 13min, with ETA
tictoc::tic()
v$distance <- calculate_shortest_route_vec(v, g, to_fuel, .progress = TRUE)
results <- v
print(results)
tictoc::toc()


# ==============================================================================
# Other Features

osmdata::available_tags("amenity")$Value
#>   [1] "animal_boarding"        "animal_breeding"        "animal_shelter"        
#>   [4] "animal_training"        "arts_centre"            "atm"                   
#>   [7] "baby_hatch"             "baking_oven"            "bank"                  
#>  [10] "bar"                    "bbq"                    "bench"                 
#>  [13] "bicycle_parking"        "bicycle_rental"         "bicycle_repair_station"
#>  [16] "bicycle_wash"           "biergarten"             "boat_rental"           
#>  [19] "boat_sharing"           "brothel"                "bureau_de_change"      
#>  [22] "bus_station"            "cafe"                   "car_rental"            
#>  [25] "car_sharing"            "car_wash"               "casino"                
#>  [28] "charging_station"       "cinema"                 "clinic"                
#>  [31] "clock"                  "college"                "community_centre"      
#>  [34] "compressed_air"         "conference_centre"      "courthouse"            
#>  [37] "crematorium"            "dancing_school"         "dentist"               
#>  [40] "dive_centre"            "doctors"                "dog_toilet"            
#>  [43] "dressing_room"          "drinking_water"         "driver_training"       
#>  [46] "driving_school"         "events_venue"           "exhibition_centre"     
#>  [49] "fast_food"              "ferry_terminal"         "fire_station"          
#>  [52] "first_aid_school"       "food_court"             "fountain"              
#>  [55] "fuel"                   "funeral_hall"           "gambling"              
#>  [58] "give_box"               "grave_yard"             "grit_bin"              
#>  [61] "hospital"               "hunting_stand"          "ice_cream"             
#>  [64] "internet_cafe"          "kindergarten"           "kitchen"               
#>  [67] "kneipp_water_cure"      "language_school"        "library"               
#>  [70] "lounge"                 "lounger"                "love_hotel"            
#>  [73] "mailroom"               "marketplace"            "monastery"             
#>  [76] "money_transfer"         "mortuary"               "motorcycle_parking"    
#>  [79] "music_school"           "music_venue"            "nightclub"             
#>  [82] "nursing_home"           "parcel_locker"          "parking"               
#>  [85] "parking_entrance"       "parking_space"          "payment_centre"        
#>  [88] "payment_terminal"       "pharmacy"               "photo_booth"           
#>  [91] "place_of_mourning"      "place_of_worship"       "planetarium"           
#>  [94] "police"                 "post_box"               "post_depot"            
#>  [97] "post_office"            "prison"                 "pub"                   
#> [100] "public_bath"            "public_bookcase"        "public_building"       
#> [103] "ranger_station"         "recycling"              "refugee_site"          
#> [106] "research_institute"     "restaurant"             "sanitary_dump_station" 
#> [109] "school"                 "shelter"                "shower"                
#> [112] "social_centre"          "social_facility"        "stage"                 
#> [115] "stripclub"              "studio"                 "surf_school"           
#> [118] "swingerclub"            "taxi"                   "telephone"             
#> [121] "theatre"                "toilets"                "townhall"              
#> [124] "toy_library"            "traffic_park"           "training"              
#> [127] "university"             "user defined"           "vehicle_inspection"    
#> [130] "vending_machine"        "veterinary"             "waste_basket"          
#> [133] "waste_disposal"         "waste_transfer_station" "water_point"           
#> [136] "watering_place"         "weighbridge"
# Created on 2024-11-07 with reprex v2.1.1