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

[Jacob van Etten]

Try to set the scaling in geoCorrection() to FALSE to get meters. (This is the default in the new version of gdistance, which will be on CRAN shortly.) See the documentation of geoCorrection() for more info.

--- On Tue, 28/6/11, Sarah Goslee <sarah.goslee at gmail.com> wrote:

> From: Sarah Goslee <sarah.goslee at gmail.com>
> Subject: Re: [R-sig-Geo] Calculating distance (km) between points by sea
> To: "Ruth Kelly" <ruthkelly123 at gmail.com>
> Cc: r-sig-geo at r-project.org
> Date: Tuesday, 28 June, 2011, 15:35
> 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
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>