Interception of routes with the number of trips

Question

I have an sf data frame with the number of trips from one bike station to another. The geometry column contains a direct route (defined by osm) between the two sides.

I want to build a map with streets colored by the gradient of the number of trips that went on.

My problem is that I have a number of trips along routes, not streets.

I use st_interception()to extract the same part from two routes and st_difference()to extract the differences.

For two lines with 10 and 15 rides here I want.

library('sf')
library('ggplot2') # dev version

route1 <- st_linestring(rbind(c(0, 0), c(1, 1), c(2, 2), c(3, 3)))
route2 <- st_linestring(rbind(c(1, 0), c(1, 1), c(2, 2), c(3, 0)))

route1 <- st_sf(id = 1, trips = 10, geometry = st_sfc(route1))
route2 <- st_sf(id = 2, trips = 15, geometry = st_sfc(route2))

# not ok as the segment (1,1 to 2,2) that is supposed to have 25 trips   only has 15 (the number of trips for the second line plotted)
ggplot(data = rbind(route1, route2)) + geom_sf(mapping = aes(color = trips)) + 
  theme(panel.grid.major = element_line(colour = 'transparent'))

# mergeRoutes gives the desired output
route <- mergeRoutes(route1, route2, init = TRUE)
ggplot(data = route) + geom_sf(mapping = aes(color = trips)) + 
  theme(panel.grid.major = element_line(colour = 'transparent'))

I wrote the mergeRoute function, it gives what I want for two routes, but cannot be extended to many, many routes.

#'
#' This function merges two routes. It returns the interscetion (if any) with the number
#' of associated count and also the remaining pars of the routes or the second route or
#' (if init) the two routes.
#'
#' @param route1 a row with id , count and geometry
#' @param route2 a row with id , count and geometry
#' @param init logical, whether to return the two routes even if there is no intersection
#'   or only the second one
#'
#' @return a data frame with 3 rows if there is an intersection, nothing otherwise.
#'   
mergeRoutes <- function(route1, route2, init = FALSE) 
{
  intersection <- st_intersection(route1$geometry, route2$geometry)
  # if the intersection is only points or is empty then the result is the two routes
  # untouched to avoid adding too many elements to the result
  if(length(intersection) != 0 & 
     !'sfc_POINT' %in% class(intersection) &
     !'sfc_MULTIPOINT' %in% class(intersection)) {
    # if intersection is a geometry with point and lines extract the lines only
    intersection <- st_collection_extract(x = intersection, type = "LINESTRING")
    count <- route1$count + route2$count
    intersection <- data.frame(id = route1$id, count = count, geometry = intersection)
    # keep the part of the initial routes that are not in the intersection
    route1_dif <- st_difference(route1$geometry, route2$geometry)
    route2_dif <- st_difference(route2$geometry, route1$geometry)
    # if one route is completely covered by the the other then it is not added to the result
    if(length(route1_dif) != 0) {
      route1 <- data.frame(id = route1$id, 
                           count = route1$count, 
                           geometry = route1_dif)
    } else {
      route1 <- NULL
    }
    if(length(route2_dif) != 0) {
      route2 <- data.frame(id = route2$id, 
                           count = route2$count, 
                           geometry = route2_dif)
    } else {
      route2 <- NULL
    }
    result <- rbind(intersection, route1, route2)
    return(result)
  } else if (init) {
    result <- rbind(route1, route2)
  } else {
    result <- route2
  }
  return(result)
}

, -, , , . , lapply() for, mac (16 RAM, 2.5 ghz), - 15 .

attemp 2000 ( ).

# To merge all the routes, each new route is compared to all the rows from the previous
# comparison. New rows are added to the resulting data frame at each step. If there is no
# intersection then the route being compared to the others is added untouched.

# initiate comparison
segment_routes <- mergeRoutes(route1 = directions %>% slice(1), 
                              route2 = directions %>% slice(2), 
                              init = TRUE)

# compute directions segmentation for all the routes
for(i in 3:nrow(directions)) {
  new_route <- directions %>% slice(i)
  # compare the new route to a the segments resulting fro mprevious comparison
  new_routes <- lapply(X = seq(nrow(segment_routes)), 
                       FUN = function(j) mergeRoutes(route1 = segment_routes %>% slice(j), 
                                                     route2 = new_route))
  new_routes <- do.call(rbind, new_routes)
  # make an sf object
  new_routes <- st_sf(new_routes, 
                      geometry = new_routes$geometry, 
                      crs = st_crs(directions))
  # add the new segemnts to the ones from the previous iteration
  segment_routes <- rbind(segment_routes, new_routes)
}

, st_intersection(), , , , , , 2 , .

sf , sp / data.table .

.

:

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.3

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_2.2.1.9000 sf_0.6-0          

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.15      class_7.3-14      withr_2.1.1.9000  plyr_1.8.4       
 [5] grid_3.4.3        gtable_0.2.0      DBI_0.7           magrittr_1.5     
 [9] e1071_1.6-8       units_0.5-1       scales_0.5.0.9000 pillar_1.2.1     
[13] rlang_0.2.0       lazyeval_0.2.1    tools_3.4.3       udunits2_0.13    
[17] munsell_0.4.3     yaml_2.1.17       compiler_3.4.3    colorspace_1.3-2 
[21] classInt_0.1-24   tibble_1.4.2     

Answer 1

, LINESTRINGS, LINESTRING - , "" "" .

data.table ( sf), .

,

1 - data.table

library(sf)
library(data.table)
library(googleway) ## for plotting. can also use ggplot2, ggmap, leaflet, mapview...

sf <- readRDS("~/Downloads/directions.rds")
sf$row_id <- 1:nrow(sf) ## for joining 

dt_routes <- as.data.table(st_coordinates(sf))

## put on the rest of the trip data
## this assumes the 'L1' value from `st_coordinates` matches the 
## `id` value from the sf_routes object
## (if not, you will need a sequential 1:nrow 'id' value to match the 
## 'L1' value)
dt_sf <- sf
st_geometry(dt_sf) <- NULL

dt_routes <- dt_routes[
  dt_sf
  , on = c(L1 = "row_id")
  , nomatch = 0
  ]

2 - 'from' 'to'

X Y, "" ""

dt_routes[, `:=`(X_to = shift(X, type = "lead"),
                 Y_to = shift(Y, type = "lead"))
          , by = L1]

3 -

dt_trips <- dt_routes[!is.na(X_to)
                      , .(n_trips = sum(count))
                      , by = .(X, Y, X_to, Y_to)]

4 - sf

/ LINESTRINGS, (.. num_trips)

dt_trips[, line_id := .I]

dt_from <- dt_trips[, .(X, Y, n_trips, line_id)]
dt_to <- dt_trips[, .(X = X_to, Y = Y_to, n_trips, line_id)]

dt_from[, line_sequence := 1]
dt_to[, line_sequence := 2]

dt_trips <- rbindlist(list(
  dt_from, dt_to
))

setorder(dt_trips, line_id, line_sequence)

## convert back to `sf` object
dt_trips <- dt_trips[, {
  geometry <- sf::st_linestring(x = matrix(c(X, Y), ncol = 2))
  geometry <- sf::st_sfc(geometry)
  geometry <- sf::st_sf(geometry)
}, by = .(line_id, n_trips)]

sf_trips <- sf::st_as_sf(dt_trips)

5 -

## applying a log-transform so the contrast shows up
sf_trips$n_trips <- log(sf_trips$n_trips)

library(googleway)
set_key("GOOGLE_MAP_KEY")

google_map(data = sf_trips) %>%
  add_polylines(stroke_colour = "n_trips", 
                stroke_opacity =1,
                stroke_weight = 3.5,
                #legend = T,
                info_window = "n_trips",
                palette = viridisLite::viridis
  )