ericpante / marmap

Import, plot and analyze bathymetric and topographic data
30 stars 7 forks source link

Issue with lc.dist(res = "dist") function #2

Closed ksamuk closed 9 years ago

ksamuk commented 9 years ago

Hi! First, this is a great package and I want to thank you for taking the time to develop it. I have run into a small issue:

I am getting strangely large values from lc.dist()

library("marmap")
library("fossil")

bat <- getNOAA.bathy(-180, 100, 30, 80, res = 10, keep = TRUE, antimeridian = FALSE)
loc <- data.frame( x = c(-149.7565, -61.9600), y = c(61.61380, 45.63197))
tr <- trans.mat(bat, min.depth = -10, max.depth = -500)
(least.cost.distance <- lc.dist(tr, loc, res="dist")[1])
[1] 3450144237

That's...a large number! Here is the lc.dist(res = "path") plotted for the same loc

lc_path_plot

So, something seemed off, so I wrote a quick function using deg.dist() from the fossil library. This function just sums the pairwise sequential distances between points in the lc.dist() path

lc_path_to_dist <- function(least.cost.path){
    path <- least.cost.path[[1]]
    dist.vec <- rep(NA, length(path[,1]))

    for (i in 2:length(path[,1])){

        long1 <- path[,1][i]
        lat1 <- path[,2][i]
        long2 <- path[,1][i-1]
        lat2 <- path[,2][i-1]

        dist.vec[i] <- deg.dist(long1, lat1, long2, lat2)

    }

    return(sum(dist.vec, na.rm = TRUE))
}

and so running:

least.cost.path <- lc.dist(tr, loc, res="path")
lc_path_to_dist(least.cost.path)
[1] 9768.091

Which seems to make more sense. I'm not exactly sure what's going on here, but I thought I'd let you know. Any input would be greatly appreciated.

besibo commented 9 years ago

Hi Kieran,

The huge value you report is due to the fact that both of your points (starting and ending points) have a depth value that lies outside the range of possible values you used to create your transition matrix. Here, you allowed a path at sea between -10m and -500m, but your starting and ending points have a positive altitude:

get.depth(bat, loc, locator = FALSE)
        lon      lat depth
1 -149.7565 61.61380    82
2  -61.9600 45.63197    35

A huge cost is thus added in the transition matrix, resulting in the seemingly incoherent value you report.

More generally, such results are observed when starting and ending points cannot be connected by staying within the given depth/altitude range specified when creating the transition matrix. This is often the case when working at large scales since points located very close to the shoreline may end up being positioned on land when using not-so-precise bathymetric data.

Here are two possible workarounds:

  1. try widening the range of depth/altitude you specify when creating the transition matrix to include the depth/altitude of both your starting and ending points. This might alter the results of the computation too strongly, and I do not encourage you to go this route.
  2. try moving slightly your starting and/or ending points so that they fall in a cell of the bathy matrix with a coherent value (here, you want a cell with a value smaller or equal to -10). You can use the get.depth() function to help you find the good spot. The function dist2isobath() can also help you find the closest locations near your starting and ending points with a suitable depth of -10m. This is probably the best solution since you should not have to move your points too far away from their "real" locations. This will thus alter the distance results only very slightly.

In your case:

new.loc <- dist2isobath(bat, loc, isobath=-10)
new.loc
  distance start.lon start.lat    end.lon  end.lat
1 60875.24 -149.7565  61.61380 -150.58333 61.23611
2 19833.29  -61.9600  45.63197  -61.76894 45.75000

The great circle distances are given in meters. So here, your distance computation would be altered by 80km at most, which seems reasonable given the scale of you map:

(least.cost.distance <- lc.dist(tr, new.loc[,4:5], res="dist")[1])
[1] 9862

I hope this helps. Please, feel free to get back to us if you can't sort it out.

Happy marmaping! Best Benoit

ksamuk commented 9 years ago

Hi Benoit,

Great, thanks for your help! That all makes sense to me. The dist2isobath() function is perfect for my application, so I think I will go with workaround 2. Thanks again!

Kieran

besibo commented 9 years ago

Great! I close this issue then. Best Benoit