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 Arctic 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)

Answer 2

To simplify the polyline you get from e.g. A* search, you can use an algorithm like Douglas-Peucker. See also this list of references: http://maven.smith.edu/~orourke/TOPP/P24.html.

Alternative idea: The usual way to apply A* would be to consider each pixel as a possible state (position), but there's no reason why you couldn't use just a subset of the pixels as possible states instead. If you make the density of states near the beginning and endpoints high, and the density of states far from either endpoint low, then you'll automatically get paths that begin and end with short, precise movements, but have long straight segments in the middle (e.g. when crossing the Pacific). If you do this, you might want to also increase the density of positions near land.

Another possible A* tweak: You can incorporate "current direction" into the state, and penalise movements that cause a change in direction. This will tend to produce long straight lines in your path. This will multiply your state space by 8, but that's probably bearable. Because you're only adding to the cost of a solution, the straight-line-to-destination heuristic you would normally use remains admissible for this new cost function, so no complications arise there.

Answer 3

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.

First create the binary sea image of the world (white: is sea, black: not sea), then erode the image. All white points after erosion are navigable. Neglecting the odd sandbank or two, of course.

As you might guess, this approach reveals a central problem in your pathfinding: Most ships have to steer quite close to land to reach a port, violating the navigation rules. However, that could be solved by starting navigation at the closest navigable sea point adjacent to a given harbour.

Answer 4

If I were you I would choose a graph theory approach. Your only problem would be to gather the edges of the database. If you have them you can use A* or Dijkstra's algo to plan the shortest route. Anyways, if I assume it right you need something similar to (Searoutefinder) right? Good luck!

Answer 5

Here is a Python approach using scgraph. In theory, you can use scgraph to solve shortest path for various network graphs using some modified Dijkstra algorithms designed to be very efficient for sparse graphs.

The package comes with some basic network graphs to get you started. There are a couple of maritime network graphs to use.

Here is an example from the docs using a slightly modified (adjusted for anti meridian support) Marnet data:

# Use a maritime network geograph
from scgraph.geographs.marnet import marnet_geograph

# Get the shortest path between 
output = marnet_geograph.get_shortest_path(
    origin_node={"latitude": 31.23,"longitude": 121.47}, 
    destination_node={"latitude": 32.08,"longitude": -81.09}
)
print('Length: ',output['length']) #=> Length:  19596.4653
print(str([[i['latitude'],i['longitude']] for i in output['coordinate_path']])) #=> [[...],...]

Depending on your graphing needs you would need to format the coordinate path accordingly for the gis system you are using.

This would produce a route like:

Note: The length produced by scgraph assumes that all non network based geograph travel has a circuity factor of 3. This is used as a way to force the optimization algorithm to choose a reasonably close network node to the provided coordinates that might not be in the network. Assuming you use an actual seaport with maritime network data, this effect should be marginal. If however you accidentally choose a port located in the middle of a landlocked area, that first leg / last leg of the journey will add 3x the great circle distance traveled.

As a side note, you might want to simplify the end route to reduce data overhead. A good option for reducing data overhead is mapshaper (web, cli).