Is there a workaround to get travel time between points (long-lat) that are not covered by googlemaps (using gmapsdistance)?

Question

I need to calculate the travel time (e.g. by car) between raster centroids (long-lat coordinates) and a reference point. However, unfortunately many of the raster centroids are not within the coverage of googlemaps and I was wondering if there is a workaround that would give me the fastest way to the reference point (by e.g. combining walking to a street and then driving to the reference point), and if so how effectively do I do this?

The problem is that I get a bias when I simply look for the next located point, take the distance to that point as walking distance and then the googlemaps driving distance from there to the reference point because the closest point might not be the one on the fastest route..

I was trying to do this using the gmapsdistance-function. Here's an example:

library(gmapsdistance)

#While this can be located and works well

> gmapsdistance(origin = "5.600451+-0.202553",
+               destination = "5.591622+-0.187677",
+               mode = "walking")
$Time
[1] 2101

$Distance
[1] 2667

$Status
[1] "OK"

#Changing the origin to other points often does not work as the points cannot be located and results in an NA-output

> gmapsdistance(origin = "7.9254948+-0.6283887",
+               destination = "5.591622+-0.187677",
+               mode = "walking")
$Time
[1] NA

$Distance
[1] NA

$Status
[1] "ROUTE_NOT_FOUND"

Many thanks!

Answer 1

One way I can think to solve this is to get a shape file of the roads in Ghana and perform geospatial operations to find the nearest road. From there you should be able to use Google's API to get the dirving distance

The steps in this solution are

1. Download shape file of roads
2. Find a 'straight' line between our origin and destination
3. Find all the roads that intersect this straight line
4. Use the nearest intersection as the new origin in the Google Directions query

---

In this example the shapefile I'm using is from here

library(sf)        ## geospatial operations
library(googleway) ## plotting and google API

---

Note:

googleway is my package, and you need a Google API key to use it

---

## Setting API keys I'll be using, one for the maps and one for directions
set_key("my_map_key"), api = "map")
set_key("my_other_api_key")

## Ghana roads
ghanaRoads <- sf::st_read("~/Downloads/gha_roads_dcw/GHA_rds_1m_dcw.shp")

## origin piont
df <- data.frame(lat = 7.9254948, lon = -0.6283887)

## destination point
dest <- data.frame(lat = 5.591622, lon = -0.187677)


google_map() %>%
    add_markers(data = df) %>%
    add_markers(data = dest) %>%
    add_polylines(ghanaRoads)

The exact approach you use here could vary. But in this example I'm using a line between the two coordinates to give us a reasonable guess at the direction of travel, and therefore where the start point should be.

## convert origin into an 'sf' object
sf_origin <- sf::st_sf(geometry = sf::st_sfc(sf::st_point(x = c(-0.6283887, 7.9254948))))

## create a line between the origin and destination
m <- matrix(c(-0.6283887, 7.9254948, -0.187677, 5.591622), ncol = 2, byrow = T)
sf_line <- sf::st_sf(geometry = sf::st_sfc(sf::st_linestring(x = m)))

## The coordinate reference system needs to match between the two for 
## spatial operations
sf::st_crs(sf_line) <- sf::st_crs(ghanaRoads)
sf::st_crs(sf_origin) <- sf::st_crs(ghanaRoads)

## find all the intersecting points
sf_intersections <- sf::st_intersection(ghanaRoads, sf_line)

google_map() %>%
    add_markers(data = df) %>%
    add_markers(data = dest) %>%
    add_polylines(data = ghanaRoads) %>%
    add_markers(data = sf_intersections)

We can use the nearest intersecting point to our origin as the starting point of our directions query.

---

Note

The sf way to find the nearest would be to use sf::st_distance, but that's dependant on lwgeom being installed on your system, but I'm having problems installing it, so I'm having to use a different method

---

I'm using a data.table function I wrote for this answer to calculate the haversine distance from each point to the origin. I then select the one with the minimum distance.

library(data.table)

coords <- matrix(unlist(sf_intersections$geometry), ncol = 2, byrow = T)

## Taking a fucntion I wrote for this answer
## https://stackoverflow.com/a/42014364/5977215
dt.haversine <- function(lat_from, lon_from, lat_to, lon_to, r = 6378137){
    radians <- pi/180
    lat_to <- lat_to * radians
    lat_from <- lat_from * radians
    lon_to <- lon_to * radians
    lon_from <- lon_from * radians
    dLat <- (lat_to - lat_from)
    dLon <- (lon_to - lon_from)
    a <- (sin(dLat/2)^2) + (cos(lat_from) * cos(lat_to)) * (sin(dLon/2)^2)
    return(2 * atan2(sqrt(a), sqrt(1 - a)) * r)
}

dt <- as.data.table(coords)
dt[, `:=`(origin_lon = df$lon, origin_lat = df$lat)]
dt[, distance := dt.haversine(origin_lat, origin_lon, V1, V2)]
## min distance
sf_nearest <- dt[order(-distance)][1, .(lon = V1, lat = V2)]
sf_nearest <- sf::st_point(c(sf_nearest$lon, sf_nearest$lat))
sf_nearest <- sf::st_sf(geometry = sf::st_sfc(sf_nearest))
sf_nearest$colour <- "green"

google_map() %>%
    add_markers(data = df) %>%
    add_markers(data = dest) %>%
    add_markers(data = sf_nearest, colour = "colour")

We can use this green marker in the directions query

orig <- sf_nearest$geometry[[1]]
orig <- as.numeric(orig)
df_orig <- data.frame(lat = orig[2], lon = orig[1])

google_map() %>%
    add_markers(df_orig)

res <- google_directions(origin = df_orig, 
            destination = dest)

## all the api results are now stored in the 'res' object.
direction_legs(res)$distance
# text  value
# 1 397 km 396829

## you can look at the route through the polygoin
df_route <- data.frame(polyline = direction_polyline(res))

google_map() %>%
    add_markers(data = df_orig) %>%
    add_markers(data = dest) %>%
    add_polylines(data = df_route, polyline = "polyline")

The dt[order(-distance)] gives us the distance from the origin to our new origin,

dt[order(-distance)][1, distance]
# [1] 1329904

This is in metres. Assume an average walking speed of 4kph and you can add it onto the overall time

---

Alternative approach

As asked for in the comments, another method for finding the nearest road is to draw a buffer around the origin and find any intersecting roads

sf_buffer <- sf::st_buffer(sf_origin, dist = 0.5)
sf::st_crs(sf_buffer) <- sf::st_crs(ghanaRoads)

google_map() %>%
  add_polylines(ghanaRoads) %>%
  add_polygons(sf_buffer)

You can then find all the lines intersecting with this buffer

sf_intersection <- sf::st_intersection(sf_buffer, ghanaRoads)

google_map() %>%
  add_markers(data = df) %>%
  add_polylines(sf_intersection)

You can use this new sf_intersection object in the distance calculation.