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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
