[R-sig-Geo] vertical profile of raster DEM along road (package raster)

[Matteo Mattiuzzi]

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 at r-project.org; omphalodes.verna at 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 at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo