Calculating distance (km) between points by sea

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Calculating distance (km) between points by sea

Ruth Kelly
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.

 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@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???

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Calculating distance (km) between points by sea

Sarah Goslee
Hi,

On Tue, Jun 28, 2011 at 9:24 AM, Ruth Kelly <[hidden email]> 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@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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]> wrote:

> From: Sarah Goslee <[hidden email]>
> Subject: Re: [R-sig-Geo] Calculating distance (km) between points by sea
> To: "Ruth Kelly" <[hidden email]>
> Cc: [hidden email]
> Date: Tuesday, 28 June, 2011, 15:35
> Hi,
>
> On Tue, Jun 28, 2011 at 9:24 AM, Ruth Kelly <[hidden email]>
> 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@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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo