How to create an interactive graph of GTFS data in R using Leaflet?

Question

I would like to create an interactive map showing the public transport routes of the city. I'm trying to do this with Leaflet in R (but am I open to alternatives, suggestions?)

Data: transport system data is in GTFS format, organized into text files (.txt), which I read in R as a data frame. *

Problem: I cannot find how to specify the identifier of each Poly line (variable shape_id) so that the graph actually matches the route of each transit line. Instead, it connects the dots in a random sequence.

Here is what I tried, so far unsuccessfully:

# Download GTFS data of the Victoria Regional Transit System tf <- tempfile() td <- tempdir() ftp.path <- "http://www.gtfs-data-exchange.com/agency/bc-transit-victoria-regional-transit-system/latest.zip" download.file(ftp.path, tf) # Read text file to a data frame zipfile <- unzip( tf , exdir = td ) shape <- read.csv(zipfile[9]) # Create base map basemap <- leaflet() %>% addTiles() # Add transit layer basemap %>% addPolylines(lng=shape$shape_pt_lon, lat=shape$shape_pt_lat, fill = FALSE, layerId =shape$shape_id) 

I would be glad to receive your comments on this subject.

* I know that you can import this data into GIS software (for example, QGIS) to create a shapefile, and then read the shapefile in R using readOGR. Robin Lovelace showed how to do it . BUT, I am looking for a clean solution R .;)

ps. Kyle Walker wrote a wonderful introduction to interactive maps in R using Leaflet . Unfortunately, it does not cover polylines in his textbook.

Answer 1

Your problem is not in the method, but in the data: note that you are loading 8 MB and that the line file that you are trying to load into the Leaflet through the brilliant one is 5 MB. As a general principle, you should always try new methods with tiny data sets first before scaling them up. This is what I do below to diagnose the problem and solve it.

Stage 1: Explore and Multiply Data

 pkgs <- c("leaflet", "shiny" # packages we'll use , "maps" # to test antiquated 'maps' data type , "maptools" # to convert 'maps' data type to Spatial* data ) lapply(pkgs, "library", character.only = TRUE) class(shape) ## [1] "data.frame" head(shape) ## shape_id shape_pt_lon shape_pt_lat shape_pt_sequence ## 1 1-39-220 -123.4194 48.49065 0 ## 2 1-39-220 -123.4195 48.49083 1 ## 3 1-39-220 -123.4195 48.49088 2 ## 4 1-39-220 -123.4196 48.49123 3 ## 5 1-39-220 -123.4197 48.49160 4 ## 6 1-39-220 -123.4196 48.49209 5 object.size(shape) / 1000000 # 5 MB!!! ## 5.538232 bytes summary(shape$shape_id) shape$shape_id <- as.character(shape$shape_id) ids <- unique(shape$shape_id) shape_orig <- shape shape <- shape[shape$shape_id == ids[1],] # subset the data 

Stage 2: Convert to Feature

data.frame this look like data.frame objects from maps?

 state.map <- map("state", plot = FALSE, fill = TRUE) str(state.map) ## List of 4 ## $ x : num [1:15599] -87.5 -87.5 -87.5 -87.5 -87.6 ... ## $ y : num [1:15599] 30.4 30.4 30.4 30.3 30.3 ... ## $ range: num [1:4] -124.7 -67 25.1 49.4 ## $ names: chr [1:63] "alabama" "arizona" "arkansas" "california" ... ## - attr(*, "class")= chr "map" 

Yes, this is similar, so we can use map2Spatial* to convert it:

 shape_map <- list(x = shape$shape_pt_lon, y = shape$shape_pt_lat) shape_lines <- map2SpatialLines(shape_map, IDs = ids[1]) plot(shape_lines) # success - this plots a single line! 

Stage 3: Join all the lines together

A for loop will make it beautiful. Please note that we use only the first 10 lines. Use 2:length(ids) for all lines:

 for(i in 2:10){ shape <- shape_orig[shape_orig$shape_id == ids[i],] shape_map <- list(x = shape$shape_pt_lon, y = shape$shape_pt_lat) shape_temp <- map2SpatialLines(shape_map, IDs = ids[i]) shape_lines <- spRbind(shape_lines, shape_temp) } 

Stage 4: The plot

Using the SpatialLines object makes the code a little shorter - in this case the first 10 lines will be built:

 leaflet() %>% addTiles() %>% addPolylines(data = shape_lines) 

Conclusion

You had to play around with the data and manipulate it before converting it to a Spatial * data type to print with the correct identifiers. maptools::map2Spatial* , unique() and a smart for loop can solve the problem.