Hi, my doubt is that extract does not sort the values according to the SpatialLine direction but sorts it somehow on its extraction order. I think it is a bug or at least a not ideal behavior in extract() on spatial lines. I think Robert can clarify this?
library(raster) # raster example in extract: r <- raster(ncol=36, nrow=18) r[] <- 1:ncell(r) cds1 <- rbind(c(-50,0), c(0,60), c(40,5), c(15,-45), c(-10,-25)) cds2 <- rbind(c(80,20), c(140,60), c(160,0), c(140,-55)) lines <- SpatialLines(list(Lines(list(Line(cds1)), "1"), Lines(list(Line(cds2)), "2") ))[1] # only 1 line! plot(r) lines(lines) x11() plot(extract(r, lines)[[1]]) # following the line the values should be sorted differently. # rasGround <- raster("rasGroundCrop.tif", crs = "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=500000 +y_0=-5000000 +ellps=bessel +towgs84=653.0,-212.0,449.0 +units=m +no_defs") Road34 <- shapefile("Profil") projection(Road34) <- projection(rasGround) plot(rasGround) plot(Road34, add = TRUE, col = "black", lwd = 3) plot(extract(rasGround, Road34)[[1]], type = "l") Matteo >>> Omphalodes Verna 04/05/13 1:24 PM >>> Dear Matteo! Thanks for suggestion. I tried your code, but it seems that it is not working. Please see my code and in attached files are files. library(sp) library(raster) library(rgdal) rasGround <- raster("rasGroundCrop.tif", crs = "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=500000 +y_0=-5000000 +ellps=bessel +towgs84=653.0,-212.0,449.0 +units=m +no_defs") Road34 <- readOGR(".", "Profil", p4s = "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=500000 +y_0=-5000000 +ellps=bessel +towgs84=653.0,-212.0,449.0 +units=m +no_defs") plot(rasGround) plot(Road34, add = TRUE, col = "black", lwd = 3) plot(extract(rasGround, Road34)[[1]], type = "l") #Not works Maybe this could be a solution but I lose direction: pointsXY <- rasterToPoints(rasGround) rasX <- rasGround rasY <- rasGround rasX[] <- pointsXY[, 1] rasY[] <- pointsXY[, 2] rasXYZ <- stack(rasGround, rasX, rasY) segment <- extract(rasXYZ, Road34)[[1]] IDX <- order(segment[, 2], segment[, 3]) segment <- segment[IDX, ] plot(segment[,1], type = "l") Thanks, OV ----- Original Message ----- From: Matteo Mattiuzzi To: r-sig-geo@r-project.org; omphalodes.ve...@yahoo.com Cc: Sent: Thursday, March 21, 2013 7:52 PM Subject: Re: [R-sig-Geo] vertical profile of raster DEM along road (package raster) Dear Omphalodes, library(raster) library(sp) set.seed(2) data(volcano) r <- raster(volcano) plot(r) l1 = cbind(c(seq(0,1 , by = 0.1)), runif(11, 0, 1)) S1 = Lines(list(Line(l1)), ID = "a") lines(S1, asp = 1, col = "red") S1<-SpatialLines(list(S1)) # I'm not practic in sp, but this at least works! extract(r,S1) Matteo >>> Omphalodes Verna 03/21/13 6:06 PM >>> Dear list. I am asking for advice for getting verticale profile of DEM along line object. I konw, it is easy for one straight line (distance and extract functions) also with using a combination of GRASS and R ( http://casoilresource.lawr.ucdavis.edu/drupal/node/375 ). Is in R possible to do this. Thanks to all. OV Here is a sample code: library(raster) library(sp) set.seed(2) data(volcano) r <- raster(volcano) plot(r) l1 = cbind(c(seq(0,1 , by = 0.1)), runif(11, 0, 1)) S1 = Lines(list(Line(l1)), ID = "a") lines(S1, asp = 1, col = "red") _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo