[R-sig-Geo] Best way to plot cross sections of discrete-valued grids

Edzer J. Pebesma e.pebesma at geo.uu.nl
Wed Apr 12 08:56:23 CEST 2006


Waichler, Scott R wrote:
> Edzer,
>
>   
>> library(sp)
>> # create 3D grid
>> xyz = expand.grid(x = 1:100, y = 1:100, z = 1:100) 
>> d = data.frame(xyz, v = sample(1:10, 1e6, replace = T))
>> gridded(d) = ~x+y+z
>> class(d)
>> summary(d)
>> # first point on line:
>> p1 = c(5,5,5)
>> # second point on line:
>> p2 = c(95,95,64)
>> rbind(p1,p2)
>> pts = sample.Line(rbind(p1,p2), 1000, "regular") # select 
>> grid elements for each of 1000 points:
>> plot(d$v[overlay(d, pts)], type = 'l')
>>
>>
>> It "misuses" sample.Line, which was written as a spsample 
>> method for Line objects; too bad that Line objects in sp are 
>> limited to exist in two dimensions. 
>>     
>
> Thank you, this example really helps me understand how to go about it.
> Here pts is along a line.  Can you suggest a way to get a set of pts
> that lie along a plane?  In my application, I want to sample and plot
> 2-D vertical cross sections ("slices") of the 3-D grid, so assuming p1
> and p2 have the same z, I need to repeat pts with different values of z.
> The only way I know how is to extract pts as a matrix and repeat using
> incremented values of z, then convert it back to a SpatialPoints object.
> There must be a more elegant way.
>
> Scott Waichler
>   
O.K., and you want to plot the different cross sections in a conditioning
levelplot? I would go like this: convert the whole thing to points (well, it
was points in my example), apply a coordinate transformation such that z
is the variable you want to cross sect, and apply

levelplot(v~x+y|z, as.data.frame(pts))

If the cross sections you want already line up with x and y, you can leave
out the coordinate transformation; the transformation makes just any
2D cross section possible. It could however be that some cross sections
don't result in nicely gridded patterns, and might need rotation before
levelplot deals with them the way you want.
--
Edzer




More information about the R-sig-Geo mailing list