New isochrone function

[rafapereirabr]

I've written a first draft of the function, available in r5r\r-package\tests_rafa\isochrone.R. As it stands, the output is a POINT sf object with all nodes the transport network and the travel time estimate between the point of interest to each node.
If we want the output to be a group of polygons, we'll need to create convex hull polygons from those points.

The funtion currently looks like this:

function

isochrone <- function(r5r_core,
                      origin,
                      destinations = NULL,
                      mode = "WALK",
                      departure_datetime = Sys.time(),
                      cutoffs = c(0, 15, 30, 45, 60),
                      max_walk_dist = Inf,
                      max_trip_duration = 120L,
                      walk_speed = 3.6,
                      bike_speed = 12,
                      max_rides = 3,
                      n_threads = Inf,
                      verbose = TRUE){

# check inputs ------------------------------------------------------------

  # check cutoffs
  checkmate::assert_numeric(cutoffs, lower = 0)



# get destinations ------------------------------------------------------------

  # if no 'destinations' are passed, use all network nodes as destination points
  if(is.null(destinations)){
    network <- street_network_to_sf(r5r_core)
    destinations = network[[1]]
  }


# estimate travel time matrix ------------------------------------------------------------

ttm <- travel_time_matrix(r5r_core=r5r_core,
                            origins = origin,
                            destinations = destinations,
                            mode = mode,
                            departure_datetime = departure_datetime,
                            max_walk_dist = max_walk_dist,
                            max_trip_duration = max_trip_duration,
                            walk_speed = 3.6,
                            bike_speed = 12,
                            max_rides = 3,
                            n_threads = Inf,
                            verbose = TRUE)

# aggregate isocrhones ------------------------------------------------------------

  # include 0 in cutoffs
  if(min(cutoffs) >0){cutoffs <- sort(c(0, cutoffs))}

  # aggregate travel-times
  ttm[, isocrhones := cut(x=travel_time, breaks=cutoffs)]


  # join ttm results to destinations
  setDT(destinations)[, index := as.character(index)]
  destinations[ttm, on=c('index' ='toId'), isocrhones := i.isocrhones]

  # back to sf
  destinations_sf <- st_as_sf(destinations)

  return(destinations_sf)
}

reprex

library(r5r)
library(ggplot2)

# build transport network
data_path <- system.file("extdata/poa", package = "r5r")
r5r_core <- setup_r5(data_path = data_path)

# load origin/point of interest
origin <- read.csv(file.path(data_path, "poa_hexgrid.csv"))[500,]

departure_datetime <- as.POSIXct("13-03-2019 14:00:00", format = "%d-%m-%Y %H:%M:%S")

# estimate travel time matrix
iso <- isochrone(r5r_core,
               origin = origin,
               mode = c("WALK", "TRANSIT"),
               departure_datetime = departure_datetime,
               cutoffs = c(0, 15, 30, 45, 60, 75, 90, 120),
               max_walk_dist = Inf,
               max_trip_duration = 120L)


ggplot() +
  geom_sf(data=iso, aes(color=isocrhones)) +
  geom_point(data=origin, color='red', aes(x=lon, y=lat)) +
  theme_minimal()

[dhersz]

I'd argue against using all network vertices as default for destinations and to force the user to pass a POLYGON/MULTIPOLYGON sf instead. Two reasons for that:

1. There are too many vertices in the network. In the relatively small sample area of PoA there are 21k vertices, while covering it with resolution 9 H3 vertices (which is pretty high resolution already) requires only 1227 hexagons.

2. We solve the problem of creating a reachable polygon from a bunch of reachable points. If the destination input is a polygon we just need to merge the reachable polygons into a single polygon and we don't need to account for error when casting points to a polygon.

[rafapereirabr]

I agree that it would be much simple if the user would pass a POLYGON/MULTIPOLYGON sf to the destinations parameter. How do other packages approach this?

[dhersz]

I've never used isochrone functions in any other packages, but apparently ropensci/opentripplanner relies in a built-in OTP isochrone API. There's no need to pass any destinations as a parameter, only origins.

[xtimbeau]

sure but you define one way or the other a resolution and you have to sample points to determine the isochrone nodes could be a good proxy for high resolution Le lun. 26 oct. 2020 à 20:47, dhersz <notifications@github.com> a écrit :

…

I've never used isochrone functions in any other packages, but apparently ropensci/opentripplanner relies in a built-in OTP isochrone API <https://github.com/ropensci/opentripplanner/blob/master/R/otp-isochrone-batch.R>. There's no need to pass any destinations as a parameter, only origins. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#123 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANA2KEJJXNRTAXT3FSM7E4DSMXG4DANCNFSM4S6V3L2A> .

[mvpsaraiva]

A few thoughts:

I like @rafapereirabr idea of using the street vertices as destinations in isochrone calculation. I was thinking on how we could implement this feature without adding an extra dependency on H3 or other grid building system, and that could be a good solution.

I also agree with @dhersz that calculations could get slow in big cities, since there are a lot of street vertices. But we would need to test it to be sure, because routing is done on the street network anyways. So, independent of how many destination points we set, all street vertices within the travel time cutoff will be visited by the algorithm. Of course, the resulting travel time matrix will be much larger for street vertices than for grid cells, so there may be an overhead in the memory exchange between Java and R.

Finally, I've checked R5 codebase and it has a built-in isochrone function that works on a grid of square cells. It also converts the reachable cells into polygons for each travel time cutoff. I just didn't have time yet to check how to make it work.

[abyrd]

Hi everyone, recently at Conveyal we were looking in more detail into R5R. It appears that as described in this issue, you are using the point-to-point routing functions to find paths one at a time to a large number of destinations. As @mvpsaraiva mentioned above, R5 has optimized methods for this use case, and in fact this is the only way Conveyal Analysis typically uses R5.

For a given origin, it computes a single shortest path tree, and then "propagates" that single tree out to a regularly spaced grid of sample points (the destinations). One such operation takes roughly the same amount of computation as a single point-to-point search. Using this one-to-many method, the number of destinations has very little effect on the computation time.

Put differently, finding paths from one origin to N destinations separately (using point-to-point routing in a loop) takes k*N units of time (roughly linear in N), but using our usual method it takes k units of time (remains roughly constant at the time required for 1 destination). This means the computations are completed thousands or millions of times faster.

Although it will require some additional work, I'd strongly encourage you to explore this option. We would have shared this earlier, but we incorrectly assumed you were using the one-to-many approach because we use it almost exclusively. In fact the point-to-point code is essentially obsolete at this point and there is a good chance we'll remove it in the near future.

[xtimbeau]

Hi Andrew, indeed, I am using massively the point to point functionality massively. What you (strongly) suggest would be a huge improvement on my way to use r5. What is the function you need to call to use this one to many function of r5 ? Thanks for your interest, Xavier Timbeau Le ven. 15 janv. 2021 à 10:44, Andrew Byrd <notifications@github.com> a écrit :

…

Hi everyone, recently at Conveyal we were looking in more detail into R5R. It appears that as described in this issue, you are using the point-to-point routing functions to find paths one at a time to a large number of destinations. As @mvpsaraiva <https://github.com/mvpsaraiva> mentioned above, R5 has optimized methods for this use case, and in fact this is the only way Conveyal Analysis typically uses R5. For a given origin, it computes a single shortest path tree, and then "propagates" that single tree out to a regularly spaced grid of sample points (the destinations). One such operation takes roughly the same amount of computation as a single point-to-point search. Using this one-to-many method, the number of destinations has very little effect on the computation time. Put differently, finding paths from one origin to N destinations separately (using point-to-point routing in a loop) takes k*N units of time (roughly linear in N), but using our usual method it takes k units of time (remains roughly constant at the time required for 1 destination). This means the computations are completed thousands or millions of times faster. Although it will require some additional work, I'd strongly encourage you to explore this option. We would have shared this earlier, but we incorrectly assumed you were using the one-to-many approach because we use it almost exclusively. In fact the point-to-point code is essentially obsolete at this point and there is a good chance we'll remove it in the near future. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#123 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANA2KENG6WHUATKCT6DSGD3S2AFBPANCNFSM4S6V3L2A> .

[abyrd]

After an email exchange with @mvpsaraiva the situation is clearer. He clarified that travel_time_matrix() is using the one-to-many approach, while the one-to-one approach is used by detailed_itineraries(). We must have been looking at the latter and jumped to conclusions. We are in the process of making more path details available in the one-to-many code. When that is finished I expect it will be possible to produce your detailed itineraries to all destinations at once, increasing the speed. That change will probably be in the next release of R5.

[mattwigway]

On the original topic of this issue, you don't want to use a convex hull - as there are closer areas that may not be accessible (i.e. isochrones may be concave). If you autogenerate a regular grid of destinations, you can use marching squares to generate vector isochrones, or with irregular points you can create a triangulated irregular network and generate contour lines from that (there may even be an R library). The isochrone drawing problem is identical to the problem of generating contour lines from elevation data, so it's well studied.

[rafapereirabr]

For the record, I've used the example in our paper to create a vignette showing how to calculate and visualize isochrones using r5r.

[rafapereirabr]

Kuan Butt has a nice post discussing Better Rendering of Isochrones from Network Graphs

[rafapereirabr]

For the record, the osrmIsochrone{osrm} function uses a really nice and simple approach. This could be a good reference.

[rafapereirabr]

Very low priority. closing this issue for now.

[rafapereirabr]

New isochrone function added to dev version in #7f50df1.

[rafapereirabr]

I have now finished a more mature version of the isochrone() function. It is basically a wrapper around the travel_time_matrix() function. It works by (1) estimating the travel times from the origin / points of interest to a set of nodes in the transport network. By default, the function considers a random sample of 80% of all nodes in the network, but users can change the sample size to improve the accuracy of the result.

Here's an example of how the function works, calculating multiple isochrones for walking from the city center of Porto Alegre and considering multiple sample sizes.

options(java.parameters = "-Xmx2G")
library(r5r)
library(ggplot2)

# build transport network
data_path <- system.file("extdata/poa", package = "r5r")
r5r_core <- setup_r5(data_path = data_path)

# load origin/point of interest
points <- read.csv(file.path(data_path, "poa_hexgrid.csv"))
origin_1 <- points[936,]

departure_datetime <- as.POSIXct(
 "13-05-2019 14:00:00",
 format = "%d-%m-%Y %H:%M:%S"
)

# estimate isochrone from origin_1
iso_0.2 <- isochrone(r5r_core,
               origin = origin_1,
               mode = c("walk"),
               departure_datetime = departure_datetime,
               cutoffs = seq(0, 100, 10),
               sample_size = .2)

iso_0.5 <- isochrone(r5r_core,
                     origin = origin_1,
                     mode = c("walk"),
                     departure_datetime = departure_datetime,
                     cutoffs = seq(0, 100, 10),
                     sample_size = .5)

iso_0.8 <- isochrone(r5r_core,
                     origin = origin_1,
                     mode = c("walk"),
                     departure_datetime = departure_datetime,
                     cutoffs = seq(0, 100, 10),
                     sample_size = .8)

iso_1.0 <- isochrone(r5r_core,
                     origin = origin_1,
                     mode = c("walk"),
                     departure_datetime = departure_datetime,
                     cutoffs = seq(0, 100, 10),
                     sample_size = 1)

iso_0.2$res <- '0.2'
iso_0.5$res <- '0.5'
iso_0.8$res <- '0.8'
iso_1.0$res <- '1.0'
iso <- rbind(iso_0.2, iso_0.5, iso_0.8, iso_1.0)



colors <- c('#ffe0a5','#ffcb69','#ffa600','#ff7c43','#f95d6a',
            '#d45087','#a05195','#665191','#2f4b7c','#003f5c')

ggplot() +
  geom_sf(data = iso, aes(fill=factor(isochrone)), color=NA) +
  scale_fill_manual(values = rev(colors) ) +
  labs(fill = 'Walking time\nin minutes') +
  facet_wrap(.~res) +
  theme_void()

[rafapereirabr]

I should add this function to the isochrone vignette in the next couple of days.

[rafapereirabr]

Done in 607ed3b.