Least-cost path = great-circle path in open ocean? (R) (gdistance)

Question

I'm using R and the gdistance package to generate the shortest routes in the ocean. In the open ocean with no obstacles, the least-cost path should equal the great-circle path. But I can't seem to get this result using gdistance (presumably because it uses a grid-based system).

Here's my approach so far (mimicking many resources on gdistance):

library(marmap) # for bathymetry data
library(gdistance)
library(raster)
library(sf)
library(leaflet)

bathy <- getNOAA.bathy(lon1 = -180, lon2 = 180, lat1 = -90, lat2 = 90, resolution = 10)

bathy_rast <- as.raster(bathy)

# Treat land as NA (non-traversable)
bathy_rast@data@values <- ifelse(bathy_rast@data@values >= 0, NA, 1) 

rast <- bathy_rast

tr_layer <- transition(rast, transitionFunction = mean, directions = 8)
tr_layer <- geoCorrection(tr_layer)

orig <- SpatialPoints(cbind(-70, 27))
dest <- SpatialPoints(cbind(-26, 33))

lc_route <- shortestPath(tr_layer, orig, dest, output = "SpatialLines")

# Inspect result
leaflet() %>%
  addTiles() %>%
  addPolylines(data = lc_route, color = "red") %>%
  addCircles(data = rbind(orig, dest))

Ideally, the line would be curved to match the great-circle path (shown above in green). Instead, it gives a somewhat unrealistic path for a ship where it goes straight, takes a turn, and then heads straight again (red).

I'm wondering if there's a way to get true shortest-distance paths in an open environment using a grid/network-based system? I still want to use a grid-based approach because I will later have routes navigate around landmasses and around obstacles (e.g., weather systems).

Answer 1

I think people have unrealistic expectations of what gdistance algorithms can do. With directions=8, any path across a free space will consist of N diagonal steps and then M horizontal (or vertical) steps. Those steps can be in any order. The path you show does the diagonal ones first. You'd get something that looks like the great circle if the diagonal and horizontal steps were interspersed.

But. The path would be the same total length, because it would still be the same number of diagonal and horizontal steps in total. It might look like a great circle, but zoom in a bit and its steps. Its like going from one corner of a city road grid to another - if you go all the way along and then all the way up the edge the distance is that same as if you wiggle along a path that is close to the diagonal (this is called the "Manhattan distance").

You can try expanding the directions parameter to 16, but that just increases the number of possible diagonal directions the path can choose. In your path case, you get a track that does a number of steps at about 22 degrees to the horizontal (ie half the slope of the first stage in the map) and then it does horizontal steps. The distance of this track will be less than the directions=8 track because it can take larger steps, but computation time is longer and its possible for the track to step over obstacles if you allow the graph of the grid to span distances.