Calculate distance raster avoiding land

Question

I'm looking for a way to calculate a raster of distances from a point, similar to what raster::distance does, but avoiding land. So, in the same way that I can make a raster of distances to a point (without any constraints) like so:

library(sp)
library(raster)


# create cooridnates sp points
col_coord <- data.frame(x      = c(-13.8309),
                        y      = c(28.9942))
coordinates(col_coord) <- ~x+y
col_coord@proj4string <- CRS(projections$WGS84)

#create base raster
baseRaster <- raster(xmn = -80, xmx =55, ymn = -70, ymx = 70, crs = CRS("+proj=longlat +ellps=WGS84"), resolution = 1)

#rasterize spatial point
colraster <- rasterize(col_coord, baseRaster)

#calculate distance
coldist <- raster::distance(colraster, doEdge = T)

#plot
plot(coldist) 

And obtain this image:

I would like to make the same raster but of "shortest distance avoiding an area" (in this case land). I'm aware of the package gdistance and I've tried to follow the advice on this post but this is as far as I've got: (Trans is a transition layer, because it takes a long time to calculate is available here.

library(gdistance)
# load transition layer (change directory if needed). This is in another CRS so we transform everything
load("data/Transition_layer.RData")

col_coord <- spTransform(col_coord, Trans@crs)
A <- accCost(Trans, col_coord)
plot(A) #this shows high cost in contintents which is good

# create a raster for water
mapWGS <- rgeos::gBuffer(rworldmap::getMap(resolution = "high"), width = 0)
water <- rasterize(mapWGS, baseRaster, mask = F, field = 1)
table(is.na(water[]))

water[water == 1] <- 0
water[is.na(water)] <- 1
water[water == 0] <- NA
plot(water)

water <- projectRaster(water, crs = Trans@crs)

# change resolution of A to match that of water
A <- resample(A, water)
A2 <- mask(A, water, inverse = F)

plot(A2) # this doesn't make that much sense. 

So, any ideas of how to get from this to a raster similar to the first one but containing distances avoiding land would be much appreciated!

Answer 1

I think you can use raster::gridDistance with an omit argument

library(raster)
xy <- data.frame(x=-13.8309, y=28.9942)
r <- raster(xmn=-120, xmx=55, ymn=-70, ymx=70, res=.25)

library(maptools)
data(wrld_simpl)

r <- rasterize(wrld_simpl, r, field= -1)
rr <- rasterize(xy, r, update=TRUE)

x <- gridDistance(rr, 1, omit=-1)
plot(x / 1000000)
points(xy)

The maps shows that it works for the Pacific, and also for the Red Sea (because wrld_simpl does not have the Panama and Suez Canals).

But because with gridDistance the path has to go through the centers of raster cells, the estimated distance is a bit longer than the actual distance.