add gtfs_isochrone fn; closes #3

DESCRIPTION

@@ -14,9 +14,12 @@ Imports:
     methods,
     Rcpp (>= 0.12.6)
 Suggests:
+    alphahull,
     hms,
     lubridate,
     knitr,
+    mapview,
+    sf,
     testthat
 LinkingTo:
     Rcpp

NAMESPACE

@@ -1,7 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(plot,gtfs_isochrone)
 export(berlin_gtfs_to_zip)
 export(extract_gtfs)
+export(gtfs_isochrone)
 export(gtfs_route)
 export(gtfs_timetable)
 importFrom(Rcpp,evalCpp)

R/RcppExports.R

@@ -17,6 +17,28 @@ rcpp_time_to_seconds <- function(times) {
     .Call(`_gtfsrouter_rcpp_time_to_seconds`, times)
 }
 
+#' rcpp_csa_isochrone
+#'
+#' Calculate isochrones using Connection Scan Algorithm for GTFS data. The
+NULL
+
+#' [deparutre_station, arrival_station, departure_time, arrival_time,
+#'     trip_id],
+#' with all entries as integer values, including times in seconds after
+#' 00:00:00. The station and trip IDs can be mapped back on to actual station
+#' IDs, but do not necessarily form a single set of unit-interval values
+#' because the timetable is first cut down to only that portion after the
+#' desired start time. These are nevertheless used as direct array indices
+#' throughout, so are all size_t objects rather than int. All indices in the
+#' timetable and transfers DataFrames, as well as startstations, are
+#' 1-based, but they are still used directly which just means that the first
+#' entries (that is, entry [0]) of station and trip vectors are never used.
+#'
+#' @noRd
+rcpp_csa_isochrone <- function(timetable, transfers, nstations, ntrips, start_stations, start_time, end_time) {
+    .Call(`_gtfsrouter_rcpp_csa_isochrone`, timetable, transfers, nstations, ntrips, start_stations, start_time, end_time)
+}
+
 #' rcpp_make_timetable
 #'
 #' Make timetable from GTFS stop_times. All stops are converted to 0-indexed

R/isochrone.R

@@ -0,0 +1,116 @@
+#' gtfs_isochrone
+#'
+#' Calculate an isochrone from a given start station, returning the list of all
+#' stations reachable to the specified `end_time`.
+#'
+#' @param gtfs A set of GTFS data returned from \link{extract_gtfs} or, for more
+#' efficient queries, pre-processed with \link{gtfs_timetable}.
+#' @param from Name of start station
+#' @param start_time Desired departure time at `from` station, either in seconds
+#' after midnight, a vector of two or three integers (hours, minutes) or (hours,
+#' minutes, seconds), an object of class \link{difftime}, \pkg{hms}, or
+#' \pkg{lubridate}.
+#' @param end_time End time to calculate isochrone
+#' @param day Day of the week on which to calculate route, either as an
+#' unambiguous string (so "tu" and "th" for Tuesday and Thursday), or a number
+#' between 1 = Sunday and 7 = Saturday. If not given, the current day will be
+#' used. (Not used if `gtfs` has already been prepared with
+#' \link{gtfs_timetable}.)
+#' @param route_pattern Using only those routes matching given pattern, for
+#' example, "^U" for routes starting with "U" (as commonly used for underground
+#' or subway routes. (Parameter not used at all if `gtfs` has already been
+#' prepared with \link{gtfs_timetable}.)
+#'
+#' @return square matrix of distances between nodes
+#' @export 
+gtfs_isochrone <- function (gtfs, from, start_time, end_time, day = NULL,
+                            route_pattern = NULL)
+{
+    if (!"timetable" %in% names (gtfs))
+        gtfs <- gtfs_timetable (gtfs, day, route_pattern)
+
+    # no visible binding note:
+    departure_time <- stop_id <- NULL
+
+    start_time <- convert_time (start_time)
+    gtfs$timetable <- gtfs$timetable [departure_time >= start_time, ]
+    if (nrow (gtfs$timetable) == 0)
+        stop ("There are no scheduled services after that time.")
+
+    stations <- NULL # no visible binding note
+    start_stns <- station_name_to_ids (from, gtfs)
+
+    stns <- rcpp_csa_isochrone (gtfs$timetable, gtfs$transfers,
+                                gtfs$n_stations, gtfs$n_trips, start_stns,
+                                start_time, end_time)
+    stns <- gtfs$stations [stns] [, stations]
+
+    stops <- gtfs$stops [match (stns, gtfs$stops [, stop_id]), ]
+    stops <- data.frame (stops [, c ("stop_name", "stop_lat", "stop_lon")])
+
+    class (stops) <- c ("gtfs_isochrone", class (stops))
+    return (stops)
+}
+
+
+#' plot.gtfs_isochrone
+#'
+#' @name plot.gtfs_ischrone
+#' @param x object to be plotted
+#' @param hull_alpha alpha value of non-convex hulls (see ?alphashape::ashape
+#' for details).
+#' @param ... ignored here
+#' @export
+plot.gtfs_isochrone <- function (x, ..., hull_alpha = 0.1)
+{
+    requireNamespace ("sf")
+    requireNamespace ("alphahull")
+    requireNamespace ("mapview")
+
+    hull <- get_ahull (x)
+
+    bdry <- sf::st_polygon (list (as.matrix (hull [, 2:3])))
+    bdry <- sf::st_sf (sf::st_sfc (bdry, crs = 4326))
+
+    x_sf <- sapply (seq (nrow (x)), function (i) {
+        sf::st_sfc (sf::st_point (as.numeric (x [i, c ("stop_lon", "stop_lat")])))
+                                })
+    x_sf <- sf::st_sf (name = x$stop_name,
+                       geometry= sf::st_sfc (x_sf, crs = 4326))
+
+    m <- mapview::mapview (x_sf, cex = 5, color = "red", col.regions = "blue",
+                           legend = FALSE)
+    mapview::addFeatures (m, bdry, color = "orange")
+}
+
+get_ahull <- function (x)
+{
+    xy <- data.frame ("x" = x$stop_lon, "y" = x$stop_lat)
+    xy <- xy [!duplicated (xy), ]
+    alpha <- 0.1
+    a <- data.frame (alphahull::ashape (xy, alpha = alpha)$edges)
+
+    xy <- rbind (data.frame (ind = a$ind1, x = a$x1, y = a$y1),
+                 data.frame (ind = a$ind2, x = a$x2, y = a$y2))
+    xy <- xy [!duplicated (xy), ]
+    xy <- xy [order (xy$ind), ]
+    inds <- data.frame (ind1 = a$ind1, ind2 = a$ind2)
+    # Wrap those indices around xy:
+    # TODO: Find a better way to do this!
+    ind_seq <- as.numeric (inds [1, ])
+    inds <- inds [-1, ]
+    while (nrow (inds) > 0)
+    {
+        j <- which (inds$ind1 == utils::tail (ind_seq, n = 1))
+        if (length (j) > 0)
+        {
+            ind_seq <- c (ind_seq, inds [j, 2])
+        } else
+        {
+            j <- which (inds$ind2 == utils::tail (ind_seq, n = 1))
+            ind_seq <- c (ind_seq, inds [j, 1])
+        }
+        inds <- inds [-j, , drop = FALSE] #nolint
+    }
+    xy [match (ind_seq, xy$ind), ]
+}

R/route.R

@@ -72,8 +72,8 @@ gtfs_route <- function (gtfs, from, to, start_time, day = NULL,
     if (nrow (route) == 0)
         stop ("No route found between the nominated stations")
 
-    route$station_name <- sapply (route$station, function (i)
-        gtfs$stops [stop_id == gtfs$stations [i] [, stations], ] [, stop_name])
+    stns <- gtfs$stations [route$station + 1] [, stations]
+    route$station_name <- gtfs$stops [match (stns, gtfs$stops [, stop_id]), ] [, stop_name]
 
     # map_one_trip maps the integer-valued stations back on to actual station
     # names. This is done seperately for each distinct trip so trip identifiers

src/RcppExports.cpp

@@ -27,6 +27,23 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// rcpp_csa_isochrone
+Rcpp::IntegerVector rcpp_csa_isochrone(Rcpp::DataFrame timetable, Rcpp::DataFrame transfers, const size_t nstations, const size_t ntrips, const std::vector <size_t> start_stations, const int start_time, const int end_time);
+RcppExport SEXP _gtfsrouter_rcpp_csa_isochrone(SEXP timetableSEXP, SEXP transfersSEXP, SEXP nstationsSEXP, SEXP ntripsSEXP, SEXP start_stationsSEXP, SEXP start_timeSEXP, SEXP end_timeSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< Rcpp::DataFrame >::type timetable(timetableSEXP);
+    Rcpp::traits::input_parameter< Rcpp::DataFrame >::type transfers(transfersSEXP);
+    Rcpp::traits::input_parameter< const size_t >::type nstations(nstationsSEXP);
+    Rcpp::traits::input_parameter< const size_t >::type ntrips(ntripsSEXP);
+    Rcpp::traits::input_parameter< const std::vector <size_t> >::type start_stations(start_stationsSEXP);
+    Rcpp::traits::input_parameter< const int >::type start_time(start_timeSEXP);
+    Rcpp::traits::input_parameter< const int >::type end_time(end_timeSEXP);
+    rcpp_result_gen = Rcpp::wrap(rcpp_csa_isochrone(timetable, transfers, nstations, ntrips, start_stations, start_time, end_time));
+    return rcpp_result_gen;
+END_RCPP
+}
 // rcpp_make_timetable
 Rcpp::List rcpp_make_timetable(Rcpp::DataFrame stop_times);
 RcppExport SEXP _gtfsrouter_rcpp_make_timetable(SEXP stop_timesSEXP) {
@@ -59,6 +76,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_gtfsrouter_rcpp_convert_time", (DL_FUNC) &_gtfsrouter_rcpp_convert_time, 1},
     {"_gtfsrouter_rcpp_time_to_seconds", (DL_FUNC) &_gtfsrouter_rcpp_time_to_seconds, 1},
+    {"_gtfsrouter_rcpp_csa_isochrone", (DL_FUNC) &_gtfsrouter_rcpp_csa_isochrone, 7},
     {"_gtfsrouter_rcpp_make_timetable", (DL_FUNC) &_gtfsrouter_rcpp_make_timetable, 1},
     {"_gtfsrouter_rcpp_csa", (DL_FUNC) &_gtfsrouter_rcpp_csa, 7},
     {NULL, NULL, 0}