Overlay spatial lines on raster

Thanks Robert,

The raster package is a magician's chest!
I get an error, though:

class(xter)

[1] "SpatialLinesDataFrame"
attr(,"package")
[1] "sp"

class(r)

[1] "RasterLayer"
attr(,"package")
[1] "raster"

lr = linesToRaster(xter, r, field = "ID_GRAFIC" )

[1] "A"
???? rows cols
[1,]??? 1?? -1
[2,]?? -1?? -1
???????? [,1]??? [,2]
[1,] 413809.3 4685685
[1]? 413786.5 4685646.9? 413810.1 4685686.4
Error in col1:col2 : NA/NaN argument
Calls: linesToRaster -> .getCols

I think because the extent of the raster is smaller than the extent of the
lines

extent(xter)

class?????? : Extent
xmin??????? : 412792
xmax??????? : 517648.1
ymin??????? : 4624794
ymax??????? : 4698657

extent(r)

class?????? : Extent
xmin??????? : 424389
xmax??????? : 512711
ymin??????? : 4636016
ymax??????? : 4685685

which I can visualize with:

plot(extent(xter))
plot(extent(r),col="red",add=T)

The problem is that I cannot find the way to clip xter to the extent of r.
pruneMap() assumes a map object and
I guess that we cannot convert from Spatial Lines to map.

I've tried

delme <- xter
delme at bbox <- bbox(r)

but still get the same problem:

lr = linesToRaster(delme, r, field = "ID_GRAFIC" )

[1] "A"
???? rows cols
[1,]??? 1?? -1
[2,]?? -1?? -1
???? [,1] [,2]
[1,]?? NA?? NA
[1]? 413786.5 4685646.9? 413810.1 4685686.4
Error in col1:col2 : NA/NaN argument
Calls: linesToRaster -> .getCols

Agus

2010/1/29 Robert J. Hijmans <r.hijmans at gmail.com>

How about something like:

library(raster)
r = linesToRaster(xter, xterpdata, field = "ID" )
v = zonal(xterpdata, r, mean)

v now shold has the mean value of the raster for each line ID and you
can merge those back


On Fri, Jan 29, 2010 at 9:54 AM, Agustin Lobo <alobolistas at gmail.com>
wrote:

Partially answering myself:
This is what I've done until now:

xter <-

readOGR(dsn="/media/Transcend/PROVI/MARICEL2/x_terRiverBasin",layer="x_ter")

class(xter)

[1] "SpatialLinesDataFrame"
attr(,"package")
[1] "sp"

xterp <- coordinates(xter)
xterp <- sapply(xterp, function(x) do.call("rbind", x))
xterp <- do.call("rbind",xterp)

as told by Roger.

Now:
require(raster)
r <- raster("test.tif")
xterpdata <- xyValues(r,xterp)
xterp2 <- data.frame(utmx=xterp[,1], utmy=xterp[,2], NO3_2007=xterpdata)
coordinates(xterp2) <- c("utmx", "utmy")

So I get to step 2, but don't see how to convert from the points
to the lines again (to get xterl). The problem is that when I did:
xterp <- sapply(xterp, function(x) do.call("rbind", x))
xterp <- do.call("rbind",xterp)

I lost the IDs of the lines. So I would need to keep the IDs of the
lines
for each point in the result of
coordinates() to be able of going back to the lines once I've put the
raster
values in the data slot of the xterp

I'll keep on, but any help would be much appreciated.

Agus

2010/1/29 Agustin Lobo <alobolistas at gmail.com>

Hi!

My goal is to colorcode a river network according to the values of
a raster map (which comes from an interpolation). The spatial lines are
imported
from a shapefile. The raster could be imported or I could run the
interpolation in R.

As what I'm planing could be rather involved, I prefer asking before:

1. Convert spatial lines to spatial points as discussed earlier this
week.
2. Overlay points on the raster and get the values for the points in a
Spatial Points Data Frame
3. Convert the Spatial Points Data Frame into Spatial Lines Data Frame
4. writeOGR() to shapefile.

The main question is: is the data frame of the points going to be
transferred to the lines
at step 3?

THnaks

Agus

? ? ? ?[[alternative HTML version deleted]]