Study

Science, Programming, Statistics & Operations Research

Work

• 2001 - Lots of (boring) things not involving transport
  
• 2009 - KiwiRail (Wellington)
  • Auckland Electrification & Track Upgrades
• 2016 - Public Transport Victoria
  • Operational Performance Analyst
• R-Ladies Talk - Visualising Spatial Data in R
• 2018 - Department of Transport
  • Transport System Analyst
  • Senior Public Transport Planner

The Problem

This project came from a meeting about access to jobs

Think about workers getting to the airport:

• Quick to drive from the West
• Free parking for staff in some cases

vs.

• Poor public transport options in the West
  • reliability
  • connections/interchanges
  • directness
  • frequency/waiting times

Modal Choice

Mode selection is a choice which a user always makes

• If driving is faster and more reliable than public
  transport for your trip, are you going to choose PT?

• Modal choice using Google Maps

• "If only we could get some data on this!"

• Famous last words: "I can do that!"

Comparing PT accessibility - Googleway

• How to compare public transport accessibility to driving accessibility at a very detailed level on the current network?

• Answer: library(googleway)!

Dave Cooley's MelbURN talk which inspired this project:
http://demos.symbolix.com.au/GooglewayRMeetupPresentation.html

Talk Plan

0. Transport/Map Jargon

1. Generating Origin Points

• Need many origin points to query

2. Running the queries in batches

• API is limited to 100 OD pairs per call

3. Extracting the results (nested lists!#@*!)

4. Using Leaflet to plot the results

5. Adjusting the scale/outlier removal

6. Where to next?

• Limitations and further work to do

Transport/Map Jargon

• PT

  • Public Transport
    (in this case including walking to stops)

• PV

  • Private Vehicle

• O-D pair

  • describes a trip between Origin and
    Destination locations

• Centroid

  • geometric (in this case) centre of a polygon

More Transport/Map Jargon

• Shapefile

  • used by GIS software to store points,
    lines and /or polygons

• Mode

  • eg driving, walking, public transport, cycling

• Mode Shift

  • moving from choosing one mode to another

• Isochrones & Choropleths - next two slides

What is an Isochrone?

What is a Choropleth?

1. Generating Origin Points

• Generate a hex grid of lots of origin points
  to query travel time to your destination from
  
  

  Need to consider:

  
  

• Excluding the River or areas with no residents
  

• Curvature of the Earth

1. Generating Origin Points

MelbourneSPDF is a spatial polygons data frame of
Metropolitan Melbourne

library(magrittr)
library(leaflet)
library(sp)
size <- 0.005
hex_points <- spsample(
  MelbourneSPDF, type = "hexagonal", cellsize = size)
length(hex_points)
[1] 36639

1. Generating Origin Points

MelbourneSPDF is a spatial polygons data frame of
Metropolitan Melbourne

library(magrittr)
library(leaflet)
library(sp)
size <- 0.005
hex_points <- spsample(
  MelbourneSPDF, type = "hexagonal", cellsize = size)
length(hex_points)
[1] 36639

1. Generating Origin Points

MelbourneSPDF is a spatial polygons data frame of
Metropolitan Melbourne

library(magrittr)
library(leaflet)
library(sp)
size <- 0.005
hex_points <- spsample(
  MelbourneSPDF, type = "hexagonal", cellsize = size)
#Want to create polygons around these points?
hex_grid <- HexPoints2SpatialPolygons(
  hex_points, dx = size)

#how do we map this in leaflet?
#call leaflet (opens interactively in Viewer window)
leaflet() %>%
  addTiles() %>% # or addTiles(url=tileURL)%>%
  addCircles(data=hex_points)

1. Generating Origin Points

Selecting the closest x points

You don't want to run up hundreds of dollars of google API costs!

How to filter the huge number of points you've generated?

• Initially I though this was a complex problem, but there's a simple solution:

  • Calculate the (rough) distance from each point to the destination
  • Sort by distance
  • Take the first x points from your data frame

• This effectively fans out origin points from the destination until
  you get to your API calls budget

• avoids messing around with selecting spatial areas to include or exclude

1. Filter the Origin Points

Calculate distance from your origin centroid & sort

#adds dist_km column to centroid_points
#from https://gis.stackexchange.com/questions/233373/distance-between-coordinates-in-r
library(data.table)
library(geosphere)
setDT(centroid_points)[, dist_km :=
  distGeo(matrix(c(lon, lat), ncol = 2),
  matrix(c(long=origin_centroid[[1]][2],lat=origin_centroid[[1]][1]), ncol = 2))/1000]
#arrange centroid_points in order of distance
centroid_points<-arrange(centroid_points,dist_km)

2. Running the queries (API Key & split)

• Set up the API Key (boring!)

  • Activate a free trial, but watch your spending: https://cloud.google.com/free

• We'll be using the Distance Matrix API, which is limited to 100 OD pairs per call

• So we need to chunk our code to send <=100 points per call

• I did this with the split() function

  • first assigning a 'tranche' number to each point

#assign a tranche number to every point
centroid_points$tranche<-factor(
  ceiling(as.numeric(rownames(centroid_points))/100)
  )
#split the list into tranches of 100
centroid_split<-split(centroid_points,centroid_points$tranche)
#result: a list of data frames

Prepare the origin points - Reformat for the API calls

library(purrr)
#fixing the format
tranches<-totalpoints/100
#initiate origins
origins<-list()
for(i in 1:tranches){
    origins[[i]]<-map2(
        (centroid_split[[i]]$lat),
        (centroid_split[[i]]$lon),
        ~c(.x,.y)
        )
}

2. Running the queries

APIcall<-function(tranche, origins, destinations, output){
#calls to google_distance by mode
  transit <- google_distance(origins = origins[[tranche]],
      destinations=destinations,       
      key = map_key,
      mode="transit",
      departure_time=depart
      )
    transit$time<-("depart.at"=paste0(depart))
    #(repeat for driving)

...

#encapsulate all output into a list chunked by tranche number
  output$transit[[tranche]] <- transit
  output$driving[[tranche]] <- driving
return(output)
}

3. Extracting the results

this is how I feel about nested lists:

(ARGHHHH!!!)     

3. Extracting the results

(Nested lists and accessors)

• You could use a really ugly for loop
  (this is definitely not what I did when I was in a hurry) ;)

• Use the googleway accessor functions:

  • this will make life slightly easier

• I would advise using functional programming

  • (do as I say, not as I do!)
    

Links:
Symbolix Googleway Accessor Vignette blog

Googleway access_result documentation

Googleway Vignette #google-distance-api

4. Using Leaflet to plot the results

library(magrittr) is your friend!

• you can build up maps layer by layer

Leaflet 101

• Background Tiles

• Adding Points Markers

• Adding Polylines

• (You can add Polygons too, and tooltips)

PTV shapefiles

Can be found at:
https://www.ptv.vic.gov.au/footer/data-and-reporting/datasets/

Or download the whole lot in one go (more than you need)
https://discover.data.vic.gov.au/dataset/public-transport-a-collection-of-ptv-datasets

Of the many files, you might be interested in these:
Paths:
Train Corridor Centreline
Metro Tram Routes
Metro Bus Routes

Stops/Stations as shapefiles:
Metro Train Stations
Metro Tram Stops
Metro Bus Stops
You can also get all the stops easily out of GTFS text files if you don't want to deal with shapefiles (see my R Ladies talk)

4. Leaflet - adding polylines (paths) and polygons

#addPolylines
leafletmap%<>% 
  addPolylines(data=metro_lines,
    weight=0.6, 
    color="#0072CE"
)
#addPolygons
leafletmap %<>% 
  addPolygons(data=SA2polygons,
    color=~pal(as.numeric(SA2_5DIG16)),
    label=~paste(SA2_NAME16,
                 SA2_5DIG16,sep=" - "),
    labelOptions = labelOptions(noHide = F), 
    fillOpacity=0.5
)

Add legend & scale bar

Legend

#where I have used this palette in leaflet
library(RColorBrewer)
pal <- colorNumeric(
    palette = "RdYlGn",reverse=TRUE,
    domain=as.numeric(outputmins$diff.mins))
leafletmap%<>% addLegend("bottomright", data=outputmins,
  pal = pal,
  values = ~diff.mins,
  title = "Minutes Difference (PT-PV)",
  opacity = 1)

Scalebar

leafletmap%<>%addScaleBar(position = "bottomleft",
    options = scaleBarOptions(maxWidth = 100,
    metric = TRUE,
    imperial = FALSE,
    updateWhenIdle = TRUE))

5. Adjusting the scale

• Points in the river
• Other outliers make it not sensible to use the default scale

• Applications: Bus Planning

• Lining up arrival not departure

• More points, closer together and for all of Melbourne

  • and multiple time samples over peak hour

• Automating querying and storing for other destinations

• Perceived reliability and number of changes

• Plotting paths taken by these predicted trips

  • googlepolylines and google directions API

• Better visualisation of the data (I'm not really a mapping person)

• See also: Conveyal and Remix