Find distance by sea from coastal point A to coastal point B

Question

I am trying to adapt the solution to this question for a slightly different purpose.

During my search for a solution, I found J. Win's nice solution here for finding a path via sea from A to B.

So I tried to adapt the code to find distances as well but did not get the expected output.

Problem 1: Scaling issue or incorrect use of function? E.g. coordinates go from north to south of Spain. The expected distance should be about 1600km~, the function outputs 66349.24

Problem 2 The required applications requires working around small islands, is there another data model that will simply slot into this when required?

library(raster)
library(gdistance)
library(maptools)
library(rgdal)
library(maps)
# make a raster of the world, low resolution for simplicity
# with all land having "NA" value
# use your own shapefile or raster if you have it
# the wrld_simpl data set is from maptools package
data(wrld_simpl)

# a typical world projection
world_crs <- crs(wrld_simpl)
world <- wrld_simpl
worldshp <- spTransform(world, world_crs)
ras <- raster(nrow=300, ncol=300)

# rasterize will set ocean to NA so we just inverse it
# and set water to "1"
# land is equal to zero because it is "NOT" NA
worldmask <- rasterize(worldshp, ras)
worldras <- is.na(worldmask)

# originally I set land to "NA"
# but that would make some ports impossible to visit
# so at 999 land (ie everything that was zero)
# becomes very expensive to cross but not "impossible" 
worldras[worldras==0] <- 999
# each cell antras now has value of zero or 999, nothing else

# create a Transition object from the raster
# this calculation took a bit of time
tr <- transition(worldras, function(x) 1/mean(x), 16)
tr = geoCorrection(tr, scl=FALSE)

# distance matrix excluding the land, must be calculated
# from a point of origin, specified in the CRS of your raster
# let's start with latlong in Black Sea as a challenging route
port_origin <- structure(c(44.206963342648436, -4.350866822197724), .Dim = 1:2)
port_origin <- project(port_origin, crs(world_crs, asText = TRUE))
points(port_origin)

# function accCost uses the transition object and point of origin
A <- accCost(tr, port_origin)

# now A still shows the expensive travel over land
# so we mask it out to display sea travel only
A <- mask(A, worldmask, inverse=TRUE)

# and calculation of a travel path, let's say to South Africa
port_destination <- structure(c(43.83115853071615, -5.134767965894603), .Dim = 1:2)
port_destination <- project(port_destination, crs(world_crs, asText = TRUE))

path <- shortestPath(tr, port_origin, port_destination, output = "SpatialLines")
t_path <-shortestPath(tr, port_origin, port_destination)
distance <-costDistance(tr, port_origin, port_destination)
# make a demonstration plot
plot(A)
points(rbind(port_origin, port_destination))
lines(path)

# we can wrap all this in a customized function
# if you two points in the projected coordinate system,
# and a raster whose cells are weighted 
# according to ease of shipping
RouteShip <- function(from_port, to_port, cost_raster, land_mask) {
  tr <- transition(cost_raster, function(x) 1/mean(x), 16)
  tr = geoCorrection(tr, scl=FALSE)
  A <- accCost(tr, from_port)
  A <- mask(A, land_mask, inverse=TRUE)
  path <- shortestPath(tr, from_port, to_port, output = "SpatialLines")
  plot(A)
  points(rbind(from_port, to_port))
  lines(path)
}

Find distance by sea from coastal point A to coastal point B

0 Answers0