[R] How to generate a smoothed surface for a three dimensional dataset?

[David Winsemius]

On Dec 5, 2013, at 4:02 PM, Jun Shen wrote:

> Thanks again, Duncan. Please allow me to ask one more question. Is it
> possible to generate a contour plot overlaying with the plot3d() plot?
> 

This is not going to give you the kewl rotation capabilities  that `rgl` offers but you might be interested in the potential of the 'plot3d' package from Karline Soetaert of NIOZ-Yerseke, The Netherlands. (She also is also part of the team that gave us the very nice ODE package, 'deSolve' and co-wrote the accompanying book. I found it interesting to read that she is a biologist, while co-authors are mathematicians.)

http://cran.rstudio.com/web/packages/plot3D/vignettes/plot3D.pdf 

See Figures 2 and 8

http://cran.rstudio.com/web/packages/plot3D/vignettes/volcano.pdf 

The second one has the catchy title: "Fifty ways to draw a volcano using package plot3D"
You may find that the examples that drew the 4 examples in Figure 7 would give you something useful.

> Jun
> 
> 
> On Thu, Dec 5, 2013 at 11:08 AM, Duncan Murdoch <murdoch.duncan at gmail.com>wrote:
> 
>> On 05/12/2013 10:33 AM, Jun Shen wrote:
>> 
>>> Hi Federico/Duncan/David/Bert,
>>> 
>>> Thank you for your thoughtful comments and it's a great learning
>>> experience. I can see the critical point here is to find a right function
>>> to make the prediction. So I was thinking to start with "loess". However
>>> the predict.loess gave me an error as follows
>>> 
>>> Error in `$<-.data.frame`(`*tmp*`, "z", value = c(0.417071766265867,
>>> 0.433916401753023,  :
>>>   replacement has 20 rows, data has 400
>>> 
>>> Here is the code I tried. Thank you for your help again!
>>> 
>>> Jun
>>> =====================================
>>> 
>>> x<-runif(20)
>>> y<-runif(20)
>>> z<-runif(20)
>>> 
>>> library(rgl)
>>> plot3d(x,y,z)
>>> 
>>> loess(z~x+y,control=loess.control(surface='direct'),
>>> span=.5,degree=2)->fit.loess
>>> 
>>> xnew <- seq(min(x), max(x), len=20)
>>> ynew <- seq(min(y), max(y), len=20)
>>> 
>>> df <- expand.grid(x = xnew, y = ynew)
>>> 
>>> df$z<-predict(fit.loess,newdata=df)
>>> 
>> 
>> After the error, use traceback() to find which function called
>> `$<-.data.frame`.  It shows that it was your final assignment
>> 
>> df$z<-predict(fit.loess,newdata=df)
>> 
>> 
>> which causes the error, because the predict function returns a matrix.  So
>> you can get the plot using
>> 
>> surface3d(xnew, ynew, predict(fit.loess,newdata=df), col="gray")
>> 
>> You may want
>> 
>> aspect3d(1,1,1)
>> 
>> afterwards; loess isn't so good at extrapolation.  Or you may want to set
>> predictions to NA outside the convex hull of your data.  (I'm not sure
>> which function is easiest to calculate that, but here's one way:
>> 
>> hullx <- x[chull(x,y)]
>> hully <- y[chull(x,y)]
>> keep <- sp::point.in.polygon(df$x, df$y, hullx, hully)
>> znew <- predict(fit.loess,newdata=df)
>> znew[!keep] <- NA
>> plot3d(x,y,z)
>> surface3d(xnew, ynew, znew, col="gray")
>> aspect3d(1,1,1)
>> 
>> 
>> 

David Winsemius
Alameda, CA, USA