1

Say I have three ships near PNG in the pacific. I want to know the shortest path between all three, while not going aground. This can be done with the ArcGIS function "cost distance" and "shortest path", but I would like to do it in R with the cool new "sf" package. Alternativly using gDistance. Anyone know how to?

Example data:

library(marmap)
library(raster)
library(sf)

#Make three ship positions near PNG
pts <- data.frame(nr=1:3, x=c(145,145,150), y=c(-10,-1,-7)) %>% 
  sf::st_as_sf(coords = c("x","y")) %>% 
  sf::st_set_crs(4326)

#Get bathymetric data from the area
papoue <- getNOAA.bathy(lon1 = 140, lon2 = 155,
                        lat1 = -13, lat2 = 0, resolution = 4)

#Turn into raster and keep only depth lower than 0 meter (water)
r <- marmap::as.raster(papoue)
plot(r)

r@data@values[r@data@values>=0] <- NA
r@data@values[r@data@values<0] <- 1

plot(r)
plot(pts, add=T)

enter image description here

1

An alternative with the gdistance package: after constructing the raster r as you do, transform your ships' positions pts into a SpatialPoints object:

pts <- sf::as(pts, "Spatial")

Then, you can use the shortestPath() function from the gdistance package. The needed input would be a raster with infinite -- NA -- cost of crossing for all landmass (which you are doing in the r raster). You need to transform the raster into a transition layer matrix with the transition() function (some extra correction might be needed, check reference):

# transforming it:
r <- transition(r, mean, directions = 8)
# Hypothetical distance from the first two ships:
distance <- shortestPath(r, pts[1,], pts[2,], output = "SpatialLines")

Good reference: http://freigeist.devmag.net/economics/683-computing-maritime-routes-in-r.html

  • Thanks, but it does not look like it is the shortest distance? imgur.com/a/EEAaC7g – Jeppe Olsen Dec 5 '18 at 13:56
  • 1
    that can be the case -- recall that diagonals are longer in distance. Either way, to make sure you account for that, you can geo-correct your transition matrix to account for that. use r <- geoCorrection(r, "c") after transforming it into a trans. matrix. The reference I gave explains this and more insightful information. – Bruno Conte Leite Dec 5 '18 at 14:49
  • Thanks. Also, it seems that using directions = 16 gives a more accurate route. – Jeppe Olsen Dec 6 '18 at 9:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.