[R-sig-Geo] Calculating distance (km) between points by sea

[Sarah Goslee]

Hi,

On Tue, Jun 28, 2011 at 9:24 AM, Ruth Kelly <ruthkelly123 at gmail.com> wrote:
> Hi,
>
> I'm trying to calculate the shortest distances by sea between points in the
> Meditteranean sea.   I've been trying to do this using a costDistance
> approach in the gdistance package.    (See code below).
>
> I have imported maps from ArcGIS in which sea has a value of 1 and land
> 10,000 there are also NA values at edges of the ascii matrix.  I have set
> the transition values to resistance so that points should choose travel by
> sea.  The projection is WCS 1984 (lat long) and the cell size (after
> aggregate command) is 0.1
>
> Hopefully the steps along the way are correct, but the output of the
> costDistance command is confusing to me.  The actually distances in km
> should be in the region of from <10 to 100km.

As a very quick first look, you're calculating distances on lat-lon
units, so you're getting distances in lat-lon units.

If I were doing it, I would reproject into m or km before calculating
the distances, since lat-lon isn't a surface distance measure. In
particular, a unit of longitude varies in length depending on the
latitude, so you really can't convert after the fact.

>  Could anybody help me to find a way of converting this?  I have tried a lot
> of variations on the code and hope it is correct.  I would appreciate
> someone looking over it for me and letting me know if it's right.
>
> Many thanks
>
> Ruth
>
> ################### code
>
> library(maptools)
>
> y <- readAsciiGrid("C:\\R\\Marine_algae\\westmed1.asc", proj4string =
> CRS("+proj=longlat + ellps=WGS84"))
>
>
> ############# vector of values from ascii grid ##########
> v1 <- y at data$westmed1.asc
>
> ########## create raster #######
>
> library(raster)
> med <- raster(y)
> setValues(med, v1)
>
> med2 <- aggregate(med, fact=10, fun=min)
>
>
> ############### g distance ##########
>
> library(gdistance)
> tr2 <- transition(med2, transitionFunction="mean", directions = 8)
> tr2 <- geoCorrection(tr2, type = "c")
> matrixValues(tr2) <- "resistance"
>
> ########## test points ##############
>
> p1 <- read.table("C:\\R\\Marine_algae\\test_points_med.csv", header = T, sep
> =",")
> p1 <- unique(p1)
> p1 <- as.matrix(cbind(p1$x, p1$y))
>
>
> cost1 <- costDistance(transition = tr2, fromCoords = p1, toCoords = p1)
>
>  cost1
>         [,1]      [,2]      [,3]      [,4]      [,5]     [,6]     [,7]
> [1,] 0.000000 2.7774957 2.7774957 2.7774957 2.0402271 2.306833 3.386609
> [2,] 2.777496 0.0000000 0.0000000 0.0000000 0.7372686 2.479186 3.576727
> [3,] 2.777496 0.0000000 0.0000000 0.0000000 0.7372686 2.479186 3.576727
> [4,] 2.777496 0.0000000 0.0000000 0.0000000 0.7372686 2.479186 3.576727
> [5,] 2.040227 0.7372686 0.7372686 0.7372686 0.0000000 1.741917 2.839459
> [6,] 2.306833 2.4791858 2.4791858 2.4791858 1.7419172 0.000000 3.106073
> [7,] 3.386609 3.5767274 3.5767274 3.5767274 2.8394588 3.106073 0.000000
>
>
> ########## ??? conversion to km???
>
>
> ___________________

-- 
Sarah Goslee
http://www.functionaldiversity.org