Thanks to the attention my paper on the cost of Somali piracy has received, a lot of people have approached me to ask how I computed the maritime routes. It is not a very difficult task using R. The key ingredient is a map of the world, that can be rasterized into a grid; all the landmass needs to be assigned an infinite cost of crossing and last but not least — one needs to compute the actual routes.
What packages do I need?
library(gdistance) library(maptools) data(wrld_simpl) library(data.table)
The package gdistance does most of the actual work of computing the routes. The wrld_simpl map provides what is needed to generate a raster.
Generating a Raster
#create a raster from shape files shp <- wrld_simpl r <- raster() r <-rasterize(shp, r, progress='text')
After the raster is generated, we can proceed by making landmass impassable for vessels.
#make all sea = -999 r[is.na(r)] <- -999 #this turns all landmass to missing r[r>-999] <- NA #assign unit cost to all grid cells in water r[r==-999] <- 1
There are a few more things to do, such as opening up the Suez Canal and some other maritime passages — one needs to find the right grid cells for this task. In the next step we can transform the raster into a transition layer matrix, that comes from the gdistance package. It is a data construct that essentially tells us how one can move from one cell to the other — you can allow diagonal moves by allowing the vessel to move into all 8 adjacent grid cells. There is also a geo-correction necessary, as the diagonals are longer distances than the straight-line moves.
tr <- transition(r, mean, directions = 8) tr <- geoCorrection(tr, "c")
Well — and thats basically it — of course, there are a few bits and pieces that need additional work — like adding heterogenuous costs as one can imagine exist due to maritime currents and so on. Furthermore, there is a whole logic surrounding the handling of the output and the storing in a local database for further use and so on.
But not to bore you with that — how can I obtain the distance between A and B? This uses Dijkstra’s Algorithm and is called through the gdistance function “shortestPath”.
AtoB <- shortestPath(tr, as.numeric(start[1:2]), as.numeric(end[1:2]), output = "SpatialLines")
Using this output, you can then generate fancy graphs such as …