Creating a trellised (faceted) thin-plate spline r

I am trying to plot a bunch of thin-plate spline response surfaces for measurements related to two continuous variables plus one discrete variable. So far, I have been subsetting the data based on the discrete variable to generate pairs of plots, but it seems to me there should be a way to create some slick trellised plots instead. It seems like this could be done by faceting heatmaps in ggplot2 with geom_tile and geom_contour, but I am stuck on

(1) how to reorganize the data (or interpret the predicted surface data) for plotting with ggplot2?

(2) syntax for creating trellised heatmaps with base graphics? or

(3) ways to use graphics from rsm to accomplish this (rsm can cope with higher-order surfaces, so I could coerce things to some extent, but plots aren't fully trellised).

Here's an example of what I've been working with so far:

library(fields)
library(ggplot2)

sumframe<-structure(list(Morph = c("LW", "LW", "LW", "LW", "LW", "LW", 
"LW", "LW", "LW", "LW", "LW", "LW", "LW", "SW", "SW", "SW", "SW", 
"SW", "SW", "SW", "SW", "SW", "SW", "SW", "SW", "SW"), xvalue = c(4, 
8, 9, 9.75, 13, 14, 16.25, 17.25, 18, 23, 27, 28, 28.75, 4, 8, 
9, 9.75, 13, 14, 16.25, 17.25, 18, 23, 27, 28, 28.75), yvalue = c(17, 
34, 12, 21.75, 29, 7, 36.25, 14.25, 24, 19, 36, 14, 23.75, 17, 
34, 12, 21.75, 29, 7, 36.25, 14.25, 24, 19, 36, 14, 23.75), zvalue = c(126.852666666667, 
182.843333333333, 147.883333333333, 214.686666666667, 234.511333333333, 
198.345333333333, 280.9275, 246.425, 245.165, 247.611764705882, 
266.068, 276.744, 283.325, 167.889, 229.044, 218.447777777778, 
207.393, 278.278, 203.167, 250.495, 329.54, 282.463, 299.825, 
286.942, 372.103, 307.068)), .Names = c("Morph", "xvalue", "yvalue", 
"zvalue"), row.names = c(NA, -26L), class = "data.frame")

sumframeLW<-subset(sumframe, Morph=="LW")
sumframeSW<-subset(sumframe, Morph="SW")

split.screen(c(1,2))
screen(n=1)
surf.teLW<-Tps(cbind(sumframeLW$xvalue, sumframeLW$yvalue), sumframeLW$zvalue, lambda=0.01)
summary(surf.teLW)
surf.te.outLW<-predict.surface(surf.teLW)
image(surf.te.outLW, col=tim.colors(128), xlim=c(0,38), ylim=c(0,38), zlim=c(100,400), lwd=5, las=1, font.lab=2, cex.lab=1.3, mgp=c(2.7,0.5,0), font.axis=1, lab=c(5,5,6), xlab=expression("X value"), ylab=expression("Y value"),main="LW plot")
contour(surf.te.outLW, lwd=2, labcex=1, add=T)
points(sumframeLW$xvalue, sumframeLW$yvalue, pch=21)
abline(a=0, b=1, lty=1, lwd=1.5)
abline(a=0, b=1.35, lty=2)

screen(n=2)
surf.teSW<-Tps(cbind(sumframeSW$xvalue, sumframeSW$yvalue), sumframeSW$zvalue, lambda=0.01)
summary(surf.teSW)
surf.te.outSW<-predict.surface(surf.teSW)
image(surf.te.outSW, col=tim.colors(128), xlim=c(0,38), ylim=c(0,38), zlim=c(100,400), lwd=5, las=1, font.lab=2, cex.lab=1.3, mgp=c(2.7,0.5,0), font.axis=1, lab=c(5,5,6), xlab=expression("X value"), ylab=expression("Y value"),main="SW plot")
contour(surf.te.outSW, lwd=2, labcex=1, add=T)
points(sumframeSW$xvalue, sumframeSW$yvalue, pch=21)
abline(a=0, b=1, lty=1, lwd=1.5)
abline(a=0, b=1.35, lty=2)

close.screen(all.screens=TRUE)

回答1:

As noted in a comment, melt() can be used to reshape the Tps() output, then it can be reformatted a bit (to remove NA's), recombined into a single data frame, and plotted. Here are plots with ggplot2 and levelplot:

library(reshape)
library(lattice)

LWsurfm<-melt(surf.te.outLW)
LWsurfm<-rename(LWsurfm, c("value"="z", "Var1"="x", "Var2"="y"))
LWsurfms<-na.omit(LWsurfm)
SWsurfms[,"Morph"]<-c("SW")

SWsurfm<-melt(surf.te.outSW)
SWsurfm<-rename(SWsurfm, c("value"="z", "X1"="x", "X2"="y"))
SWsurfms<-na.omit(SWsurfm)
LWsurfms[,"Morph"]<-c("LW")

LWSWsurf<-rbind(LWsurfms, SWsurfms)

LWSWp<-ggplot(LWSWsurf, aes(x,y,z=z))+facet_wrap(~Morph)
LWSWp<-LWSWp+geom_tile(aes(fill=z))+stat_contour()
LWSWp

or: levelplot(z~x*y|Morph, data=LWSWsurf, contour=TRUE)

回答2:

require(rgl)
open3d()
plot3d
surface3d(surf.te.outSW$x, surf.te.outSW$y, surf.te.outSW$z, col="red")
surface3d(surf.te.outLW$x, surf.te.outLW$y, surf.te.outLW$z, col="blue")
decorate3d()
      rgl.snapshot("OutRGL.png")

Another version where I scaled the x and y values by a factor of 10 and rotated to "look through" the gap. If this were your choice you might want to look at ?scaleMatrix