How to generate a smoothed surface for a three dimensional dataset?

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

Jun


On Thu, Dec 5, 2013 at 11:08 AM, Duncan Murdoch
<murdoch.duncan at gmail.com <mailto: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)