algorithm - Find path by sea from coastal point A to coastal point B

Question

I have the seemingly tricky challenge of trying to work out a path, by sea, from one sea port to another sea port. The eventual aim is to plot this on a Google (or Bing) Map as a polyline.

The path needs to:

• Be plausible, as a ship can't go over land (obviously)
• Not run too close to the coast line. Ships can't go that far near to shore
• Not be too complex. It's going to be plotted on Google Maps so a 2000 point polyline wont do.
• Be the shortest, but not at the expense of the above three points

So, my first idea was obtain data on coast lines around the world. Such a thing is available here. Unfortunately it's incomplete however. OpenStreetMap shows this data and shorelines for things like Caribbean islands are missing.

I also thought about Geocoding (not reliable enough plus I would burn through thousands of requests trying to plot a route)

My next idea was to somehow use Google Maps and test if a point is blue or not. GMaps.NET, a great .NET Mapping component, allowed me to achieve this by creating a bitmap of what it renders and testing the color of a pixel.

First problem is that the accuracy of this hit testing is only as good as the resolution image of the image I test against. For ports close to each other, this is fine for ports further away, the accuracy suffers.

Second problem, assuming I use some kind of 'blue pixel testing' method, is what algorithm is right for finding a route. The A* algorithm looks promising, but I'm not sure how to push the path 'out' from the being to near to the coast. Nor how to reduce complexity of the polyline.

So ... any input: ideas, thoughts, links, sample code, etc would be welcome. Thanks.

(I should add that this is for a travel site. Accuracy isn't too important, I'm not directing shipping or anything)

Answer 1

Here's a solution using R. Could be much refined, simply by making certain regions expensive or cheap for the shortest path algorithm. For example exclude the Artic Ocean, allow major rivers/canals, make known shipping routes preferred for travel.

library(raster)
library(gdistance)
library(maptools)
library(rgdal)

# 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(30, 40), .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(20, -35), .Dim = 1:2)
port_destination <- project(port_destination, crs(world_crs, asText = TRUE))

path <- shortestPath(tr, port_origin, port_destination, output = "SpatialLines")

# 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)
}

RouteShip(port_origin, port_destination, worldras, worldmask)