convert x,y,z data into a grid

Dear list members,
I'm looking for a function in R like the function r.in.xyz in GRASS

6.2). I have x, y, z data frame (UTM) with z equals depth and would 
like to make a 3dperspective plot with rgl. My dataframe has 2.5 
milion rows, so I think it's necessary to convert the data to a grid 
with something like mean depth per grid cell. I know I've seen an 
answer to a similar question somewhere on one of the lists, but I just

don't seem to be able to find it again. Thanks for any help.
Erik Meesters
 

Dear Erik,

In my opinion you beter first create an interpolated grid (e.g. with
kriging). Then plot this interpolated grid. That will reduce the number
of points dramatically. A 100 by 100 grid will probably do.

HTH,

Thierry


------------------------------------------------------------------------
----
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
methodology and quality assurance
Gaverstraat 4
9500 Geraardsbergen
Belgium
tel. + 32 54/436 185
Thierry.Onkelinx at inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to
say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of
data.
~ John Tukey

-----Oorspronkelijk bericht-----
Van: r-sig-geo-bounces at stat.math.ethz.ch
[mailto:r-sig-geo-bounces at stat.math.ethz.ch] Namens Meesters, Erik
Verzonden: woensdag 7 mei 2008 14:23
Aan: r-sig-geo at stat.math.ethz.ch
Onderwerp: Re: [R-sig-Geo] convert x,y,z data into a grid

Michal,
thanks for the suggestions. If I try the wireframe I get a warning
message

"Warning message:
In nx * ny : NAs produced by integer overflow"

I haven't been able to figure out why that is, but could have something
to do with the 2.5mil * 2.5mil number of rows....

My data are not in a grid format, so the second option will not work.
The data are simply UTM positions and depth values from a cruise
covering approximately an area of 10 square kilometers. In Grass there
is a function that reads in the dataframe and converts the data to a
grid which I think is easier to work with in many of the different
plotting routines.

Hope anyone can help.
Erik



-----Original Message-----
From: Michal Gallay [mailto:mgallay01 at qub.ac.uk]
Sent: Wednesday, May 07, 2008 11:48
To: Meesters, Erik
Subject: Re: [R-sig-Geo] convert x,y,z data into a grid

On May 7 2008, Meesters, Erik wrote:

Dear list members,
I'm looking for a function in R like the function r.in.xyz in GRASS

(vs.

6.2). I have x, y, z data frame (UTM) with z equals depth and would
like to make a 3dperspective plot with rgl. My dataframe has 2.5
milion rows, so I think it's necessary to convert the data to a grid
with something like mean depth per grid cell. I know I've seen an
answer to a similar question somewhere on one of the lists, but I just

don't seem to be able to find it again. Thanks for any help.
Erik Meesters

Hello Erik,

you could try something like the following code with lattice package.
You don't need to convert the dat to grid for this.

require(lattice)

mydata <- read.table(...)

wireframe(z ~ x * y, data = mydata,
          scales = list(arrows = FALSE),
          drape = TRUE, colorkey = TRUE,
          screen = list(z = 30, x = -60))

or with simple graphics function 'persp' but for this you need to
convert the data to matrix/grid. The following code presumes your x,y,z
coordinates are on a grid. It will not work if they are not or at least
I am not aware of any other approach.

require(sp)

coordinates(mydata) <- ~ x+y
gridded(mydata) <- TRUE
m <- as.matrix(mydata)
exag <- 3
spacing <- 5
theta_val <- 120
phi_val <- 35
ltheta_val <- 20
lphi_val <- 10
z <- exag * m      # Exaggerate the relief
x <- spacing * (1:dim(m)[1])
y <- spacing * (1:dim(m)[2])

par(c(0.5,0.5,0.5,0.5))
persp(x, y, z, theta = theta_val, phi = phi_val, col = "gray", scale =
FALSE, ltheta = ltheta_val, lphi=lphi_val, shade = 0.75, border = NA,
box = TRUE, r=600)

Both of them worked for my 2mil points data.

I hope this helps.

Michal

--
Michal Gallay

Postgraduate Research Student
School of Geography, Archaeology and Palaeoecology Queen's University
Belfast BT7 1NN Northern Ireland

Tel: +44(0)2890 273929
Fax: +44(0)2890 973212
email: mgallay01 at qub.ac.uk
www: www.qub.ac.uk/geog