R/make_isochrone.R

In Chrisjb/RDistanceMatrix: Useful Functions For Building Isochrones And Related Analysis

#' Make Time Isochrone
#'
#' Generates a time isochrone for travel times to/from a specified location.
#'
#' @param site either a data.frame with the columns 'lat' and 'lng' or a character string to be geocoded.
#' @param direction either 'in' or 'out'.'in' generates an isochrone from multiple origins to the site. 'out' goes from the site to multiple destinations.
#' @param time the travel time in minutes for the isochrone extents
#' @param detail parameter used for method=google and takes 'low', 'medium' or 'high' level of detail in the isochrone. High will produce the most granular detail but will use more API credits.
#' @param multiplier used to calculate the extents of the google isochrone. It should be tuned using method='google_guess'. In areas with a low speed of traffic, lower multiplier values will be appropriate (~0.7). When the mode is set to walking, a recommended parameter is ~0.1. The interpretation of this value is that a multiplier of 1 means that the maximum distance able to be travelled in 1 minute is 1km.
#' @param mode a character string for the mode of travel. Possible values are 'driving', 'cycling', 'transit', or 'walking'
#' @param method eithing 'google' or 'mapbox' method for generating isochrones. Read the vingette for more details. 'google_guess' is also available to tune the multiplier parameter.
#' @param departing  departure time can be set if we know the departure time and want to take into account traffic in the network for driving or transit modes.
#'     Departing should always be set to FALSE for walking or cycling modes. If required, departure time should be specified in the format 'YYYY-MM-DD HH:MM:SS', in the local time of the origin country
#' @param tz if departure time is set, timezone defaults to the user's timezone as per `Sys.timezone`. Where the timezone of the `site` differs from the user timezone, specify the `site` timezone here.
#'     tz is a character string that must be a time zone that is recognized by the user's OS.
#' @param init_grid_size used for method='google' - adjust if we want to pick up more points in the 'google_guess' stage before adding detail. Useful for longer travel times (>60mins).
#' @param google_api_key the google maps API key. This can be generated from the google cloud console and set using 'set_google_api'
#' @param mapbox_api_key the mapbox maps API key. This can be generated by signing up to a mapbox account and set using 'set_mapbox_api'
#'
#' @return an sf polygon
#'
#' @import sf lubridate leaflet
#' @importFrom dplyr mutate select rename if_else case_when bind_rows
#' @importFrom maptools ContourLines2SLDF
#' @importFrom utils URLencode
#' @importFrom httr GET content http_error http_status
#'
#' @examples
#' transit_isochrone <- make_isochrone(site = 'London Bridge Station, London, UK', time = 30,
#'   direction ='out',  multiplier = 0.7, mode = 'transit')
#' driving_isochrone_mapbox <- make_isochrone(site = 'London Bridge Station, London, UK',
#'  time = 30, method = 'mapbox', mode = 'driving')
#'
#' library(leaflet)
#' leaflet(driving_isochrone_mapbox) %>%
#'   addTiles() %>%
#'   addPolygons()
#'
#' @export
make_isochrone <- function(site, time, direction = c('out','in'),
                           detail = 'medium', multiplier = 0.8,
                           mode= c('driving','walking', 'cycling', 'transit'),
                           departing = F, method = c('google', 'mapbox', 'google_guess'),
                           tz = Sys.timezone(),
                           init_grid_size = 5,
                           google_api_key = Sys.getenv('google_api_key'),
                           mapbox_api_key = Sys.getenv('mapbox_api_key')){

  # handle match.arg
  mode <- match.arg(mode)
  method <- match.arg(method)
  direction <- match.arg(direction)
  args <- as.list(match.call())[-1]
  form <- formals()
  form[names(args)] <- args
  form['mode'] <- mode
  form['method'] <- method
  form['direction'] <- direction

  # validity checks...
  do.call(check_request_valid, form)

  # geocode site if lat lng not already given
  if(is.character(site)) {
    message('Geocoding: "', site,'" if you entered precise co-ordinates, please specify site as a data frame containing the columns "lat" and "lng"')
    if(method %in% c('google', 'google_guess')){
      site <- geocode(site, api_key = api_key)[1,]
    } else {
      site <- geocode_mapbox(site, api_key = api_key)[1,]
    }
  }

  # make isochrone
  if(method == 'mapbox'){
    coordinates <- paste0(site$lng,',',site$lat)
    if(mode == 'transit'){
      stop("transit isochrones are not available with method mapbox")
    }
    mapbox_url <- glue::glue('https://api.mapbox.com/isochrone/v1/mapbox/{mode}/{coordinates}?contours_minutes={time}&polygons=true&access_token={api_key}')
    mapbox_url_secret <- gsub(api_key, "SECRET", mapbox_url)
    message(glue::glue('fetching isochrone from url: {mapbox_url_secret}'))
    mb_res <- httr::GET(mapbox_url)
    if(httr::http_error(mapbox_url)){
      stop(httr::content(mb_res)$message)
    }
    poly <- sf::read_sf(mapbox_url)
    if(all(sf::st_is_valid(poly))==F){
      poly <- sf::st_make_valid(poly)
    }
    return(poly)
  }

  message('drawing initial isochrone...')
  distance_guess <- time * multiplier
  extents_guess <- approx_grid(site$lat, site$lng, distance_guess)

  # initial grid size of 5x5, which will make the max dimension requested be 25
  lats_init <- seq(from= -extents_guess$lat, to = extents_guess$lat, length.out = init_grid_size) + site$lat
  lngs_init <- seq(from= -extents_guess$lng, to = extents_guess$lng, length.out = init_grid_size) + site$lng

  df_init <- expand.grid(lats_init, lngs_init)
  names(df_init) <- c('lat', 'lng')

  df_init <- dplyr::mutate(df_init, latlng=paste0(lat,',',lng))

  # set up for direction = out
  if(direction == 'out') {
    dists_init <- distance_to_destinations(origin = paste0(site$lat,',',site$lng), dest = df_init$latlng, mode = mode,
                                      departing = departing, tz = tz,  api_key = api_key)
    # set up for direction = in
  } else if(direction == 'in'){
    dists_init <- distance_from_origins(origin = df_init$latlng, dest = paste0(site$lat,',',site$lng),
                                   mode = mode, departing = departing, tz = tz, api_key = api_key)
  }

  df_init <- dplyr::mutate(df_init, minutes = as.numeric(dists_init$time)/60)


  # make distance matrix
  m_init <- matrix(df_init$minutes, nrow = sqrt(nrow(df_init)),
              dimnames = list(unique(df_init$lat),
                              unique(df_init$lng)),
              byrow = T)

  # if missing value, make it arbitrarily high
  m_init[is.na(m_init)] <- 2*time

  tryCatch({
    contour_init <- contourLines(x=unique(df_init$lng), y= unique(df_init$lat), z=m_init, levels = seq(from = 0, to = time, by= time/5 ))
  },
  error =function(e){
    message(glue::glue("Unable to find an isochrone. It's likely the multiplier is set too high. Original error message: {e}"))
  })


  shp_init <- maptools::ContourLines2SLDF(contour_init)

  shp_init <- shp_init %>% sf::st_as_sf() %>%
    sf::st_set_crs(4326) %>% sf::st_cast('POLYGON') %>%
    sf::st_make_valid() %>% sf::st_union()


  if(method == 'google_guess'){
    chrone_map <- leaflet::leaflet(shp_init) %>%
      leaflet::addPolygons(col='blue',fillColor = 'yellow') %>%
      leaflet::addTiles() %>%
      leaflet::addMarkers(data=site) %>%
      leaflet::addCircles(data=df_init,lng = ~lng,lat=~lat,col='grey',radius = 0.2,popup=~paste0(lng,', ',lat,'  TIME:',minutes))
    credits <- dplyr::if_else(departing == F, length(df_init$latlng)*0.005, length(df_init$latlng)*0.01)
    message(glue::glue('Google API elements used: {length(df_init$latlng)} (\u00A3{credits} credits)'))
    return(chrone_map)
  }

  # improve detail on initial shape
  message('adding detail to initial isochrone...')

  # get bounding box of our google maps initial poly to calibrate the buffers
  bbox <- sf::st_bbox(shp_init)
  left <- data.frame(lng = bbox[['xmin']], lat = bbox[['ymin']]) %>%
    sf::st_as_sf(coords = c('lng', 'lat'), crs = 4326)
  right <- data.frame(lng = bbox[['xmax']], lat = bbox[['ymin']]) %>%
    sf::st_as_sf(coords = c('lng', 'lat'), crs = 4326)
  top <- data.frame(lng = bbox[['xmin']], lat = bbox[['ymax']]) %>%
    sf::st_as_sf(coords = c('lng', 'lat'), crs = 4326)
  bottom <- data.frame(lng = bbox[['xmin']], lat = bbox[['ymin']]) %>%
    sf::st_as_sf(coords = c('lng', 'lat'), crs = 4326)
  width <- as.numeric(st_distance(left, right))
  height <- as.numeric(st_distance(top, bottom))

  # set buffer to 15% of max(height,width)
  buffer_dist <- max(height, width) * 0.15


  # find key buffer region where shp_range is the possible region for isochrone extents
  shp_max <- shp_init %>%
    sf::st_transform(27700) %>%
    sf::st_buffer(buffer_dist) %>%
    sf::st_transform(4326)

  shp_min <- shp_init %>%
    sf::st_transform(27700) %>%
    sf::st_buffer(-buffer_dist) %>%
    sf::st_transform(4326)

  shp_range <- suppressMessages(
    shp_max %>% sf::st_difference(shp_min)
  )


  # create grid of points over our bounding box
  grid_bbox <- sf::st_bbox(shp_max)

  pts <- dplyr::case_when(detail == 'min' ~ 15,
                          detail == 'low' ~ max(floor(max(height, width)/2500),15),
                          detail %in% c('med','medium') ~ max(floor(max(height, width)/1000), 30),
                          detail == 'high' ~ max(floor(max(height, width)/500),40))

  tolerence <- max(height, width)/pts

  lats <- seq(from=grid_bbox$ymin, to = grid_bbox$ymax, length.out = pts)
  lngs <- seq(from=grid_bbox$xmin, to = grid_bbox$xmax, length.out = pts)
  grid <-  expand.grid(lat = lats, lng = lngs)

  grid <- grid %>%
    sf::st_as_sf(coords = c('lng','lat'),crs= 4326, remove = F)


  # intersect grid with possible area range
  sf::st_agr(grid) ='constant'
  donut <- suppressMessages(
    grid %>%
    sf::st_intersection(shp_range)
  )

  # create matrix from our grid of points
  pts_mat <- grid  %>%
    dplyr::mutate(latlng = paste0(lng,',',lat))

  # create points inside donut hole and set travel time to just less than the max time
  inner_pts <- suppressMessages(
    grid %>%
    sf::st_intersection(shp_min) %>%
    dplyr::mutate(latlng = paste0(round(lat,5),',',round(lng,5)),
           minutes = time -1)
    )

  # create points in donut ring for geocoding
  df <- donut %>%
    dplyr::mutate(latlng = paste0(lat,',',lng))

  elements <- length(df_init$latlng)+length(df$latlng)
  credits <- dplyr::if_else(departing == F, elements*0.005, elements*0.01)

  if(credits >= 10) {
    message(glue::glue("Large Query Warning: request will use {elements} elements (\u00A3{credits} credits)."))
    response <- readline(prompt=paste0("Continue? (y/n)"))

    if(!response %in% c('y', 'Y')) {
      stop('Use a lower detail setting')
    }

  }
  # set up for direction = out
  if(direction == 'out') {
    dists <- distance_to_destinations(origin = paste0(site$lat,',',site$lng), dest = df$latlng, mode = mode,
                                           departing = departing, tz = tz, api_key = api_key)
    # set up for direction = in
  } else if(direction == 'in'){
    dists <- distance_from_origins(origin = df$latlng, dest = paste0(site$lat,',',site$lng),
                                        mode = mode, departing = departing, tz =tz,  api_key = api_key)
  }

  df <- dplyr::mutate(df, minutes = as.numeric(dists$time)/60)

  # join matrices
  distance_df <- suppressMessages(
    pts_mat %>%
    sf::st_join(df) %>%
    dplyr::select(lng=lng.x, lat=lat.x, latlng=latlng.x, minutes) %>%
    sf::st_join(inner_pts,left = T) %>%
    dplyr::mutate(minutes = dplyr::if_else(!is.na(`minutes.x`), `minutes.x`, `minutes.y`)) %>%
    dplyr::select(lng = lng.x, lat = lat.x, latlng = latlng.x, minutes)
    )

  # make distance matrix
  m <- matrix(distance_df$minutes, nrow = sqrt(nrow(distance_df)),
              dimnames = list(unique(distance_df$lat),
                              unique(distance_df$lng)),
              byrow = T)


  # if missing value, make it arbitrarily high
  m[is.na(m)] <- time*1.1

  contour <- contourLines(x=unique(distance_df$lng), y= unique(distance_df$lat),
                       z=m, levels = seq(from = 0, to = time, by = time/5 ))

  shp <- maptools::ContourLines2SLDF(contour)

  shp <- shp %>% sf::st_as_sf() %>% sf::st_set_crs(4326) %>%
    sf::st_cast('POLYGON') %>% sf::st_make_valid() %>%
    sf::st_union() %>%
    sf::st_sf()

  message(glue::glue('Google API elements used: {elements} (\u00A3{credits} credits). Isochrone generated to accuracy of {round(tolerence)}m'))
  if(all(sf::st_is_valid(shp))==F){
    shp <- sf::st_make_valid(shp)
  }
  return(shp)

}