[R-sig-Geo] Overlay spatial lines on raster

[Robert J. Hijmans]

Hi Agus, thank you. You were right about the cause of the problem ---
I patched it in version 0.9.9-5, now on R-Forge (and should be
available for automatic install within 24 hours); the function could
use some optimization for speed but at least it works now for
SpatialLines* with a larger extent then the target RasterLayer.

# simple example
library(raster)
library(maptools)
data(wrld_simpl)
r = raster(xmn=-10, xmx=20, ymn=37, ymx=43, ncol=360, nrow=72)
r = linesToRaster(wrld_simpl, r, progress='text')
plot(r)
plot(wrld_simpl, add=T)


# example using zonal as I proposed in your original question, but
here with polygonsToRaster
# What is (roughly) the mean latitude of each country?

x = raster()
# res(x) = 0.1    # higher res, slower, but more small countries are included
x = polygonsToRaster(wrld_simpl, x, progress='text')
y = raster(x)
y[] = rep(yFromRow(y, 1:nrow(y)), each=ncol(y))
z = zonal(y, x, mean)
res= data.frame(wrld_simpl[z[,'zone'], 'NAME'], lat=z[,2])
res[order(res[,2]), ]

Robert



On Sat, Jan 30, 2010 at 5:03 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
> 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]]
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo at stat.math.ethz.ch
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>
>