[R-sig-Geo] convert x,y,z data into a grid

Wed May 7 16:08:18 CEST 2008

This would be my way, with kernel smoothing from spatstat function smooth.ppp:

Let's say that "cosa" is a dataframe with three columns x, y and z .

library(rgl)

library(ecespa) # for function "haz.ppp"; it requires (and load)  spatstat

cosa.ppp <- haz.ppp (cosa) # makes ppp object 
"easily". You may consider also functions "ppp" or "as.ppp" in spatstat.

cosa.s <- smooth.ppp (cosa.ppp)  #  spatial 
smoothing of numeric values observed

cosa.p <- persp.im (cosa.s) # 3D transformation 
matrix returned by persp.default

persp3d(cosa.p) # your perspective plot

Regards, Marcelino

At 14:33 07/05/2008, ONKELINX, Thierry wrote:
>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
>
>
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>--
>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
>
>_______________________________________________
>R-sig-Geo mailing list
>R-sig-Geo at stat.math.ethz.ch
>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>_______________________________________________
>R-sig-Geo mailing list
>R-sig-Geo at stat.math.ethz.ch
>https://stat.ethz.ch/mailman/listinfo/r-sig-geo

________________________________

Marcelino de la Cruz Rot

Departamento de  Biología Vegetal
E.U.T.I. Agrícola
Universidad Politécnica de Madrid
28040-Madrid
Tel.: 91 336 54 35
Fax: 91 336 56 56
marcelino.delacruz at upm.es