Hi, On Tue, Jun 28, 2011 at 9:24 AM, Ruth Kelly <ruthkelly...@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@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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo