display 3d plots

x=seq(-pi,pi,len=50)
y=x;
f=outer(x,y,function(x,y)cos(y)/(1+x^2));
f[1:5, 1:5]

# The contour() function produces a contour map and is used for three
# dimensional data (like mountains). You feed it three inputs. The first
# is a vector of the x values, the second a vector of the y values and
# the third is a matrix with each element corresponding to the Z value
# (the third dimension) for each pair of (x,y) coordinates. Just like
# plot there are many other things you can feed it to. See the help file.

contour(x,y,f)
contour(x,y,f,nlevels=15)

# persp() works the same as image() and contour() but it actually
# produces a 3d plot. theta an phi control the various angles you can look
# at the plot from.
persp(x,y,f)
persp(x,y,f,theta=30)
persp(x,y,f,theta=30,phi=20)



########################################################################
## another example of 3d plot from my personal reserach, use rgl library
########################################################################
# 3D visualization device system

library(rgl);
data(volcano)
dim(volcano)

peak.height <- volcano;
ppm.index <- (1:nrow(volcano));
sample.index <- (1:ncol(volcano));

zlim <- range(peak.height)
zlen <- zlim[2] - zlim[1] + 1
colorlut <- terrain.colors(zlen) # height color lookup table
col <- colorlut[(peak.height-zlim[1]+1)] # assign colors to heights for each point
open3d()

ppm.index1 <- ppm.index*zlim[2]/max(ppm.index);
sample.index1 <- sample.index*zlim[2]/max(sample.index)

title.name <- paste("plot3d ", "volcano", sep = "");
surface3d(ppm.index1, sample.index1, peak.height, color=col, back="lines", main = title.name);
grid3d(c("x", "y+", "z"), n =20)

sample.name <- paste("col.", 1:ncol(volcano), sep="");
sample.label <- as.integer(seq(1, length(sample.name), length = 5));

axis3d('y+',at = sample.index1[sample.label], sample.name[sample.label], cex = 0.3);
axis3d('y',at = sample.index1[sample.label], sample.name[sample.label], cex = 0.3)
axis3d('z',pos=c(0, 0, NA))

ppm.label <- as.integer(seq(1, length(ppm.index), length = 10));
axes3d('x', at=c(ppm.index1[ppm.label], 0, 0), abs(round(ppm.index[ppm.label], 2)), cex = 0.3);

title3d(main = title.name, sub = "test", xlab = "ppm", ylab = "samples", zlab = "peak")
rgl.bringtotop();

special symbols on R plot

Sometimes we need put special symbols, like Greek letters and italic type, on the plots. Sometimes, we also want to put mathematical annotation on the plot. In R we can use the function expression() do this job:

Sample codes

xlab.name = expression(paste(italic(vti), Delta, sep=""))
ylab.name = expression(mu * "ml")
main.name = expression(paste(plain(sin) * phi))
plot(0, 0, xlab=xlab.name, ylab=ylab.name, main=main.name, xlim=c(-pi, pi), ylim=c(-1.5, 1.5), axes=FALSE)
axis(1, at = c(-pi, -pi/2, 0, pi/2, pi), labels = expression(-pi, -pi/2, 0, pi/2, pi))
axis(2)
box()
text(-pi/2, 0, expression(hat(alpha) == (X^t * X)^{-1} * X^t * y))
text(pi/2, 0, expression(paste(frac(1, sigma*sqrt(2*pi)), plain(e)^{frac(-(x-mu)^2, 2*sigma^2)}, sep="")), cex = 1.2)

### #######################
PS: for the subscript expression,
We should use
> ylab.name = expression(sigma[21])

############################
There are webpages talking about this issue in details. Mathematical Annotation in R: http://rweb.stat.umn.edu/R/library/grDevices/html/plotmath.html
##############################
others:
main=expression(paste(italic(Sch9), Delta, " (18s and 5.8s)")),

However, if we want to use a variable in the expression, we should use substitute:
e.g.

n <- 20
plot(0, 0, main = substitute(paste(n[i], " = ", k), list(k = n)));
i=2;
range.name = substitute(paste(italic(Sch9), Delta, " ", p, "/", k, " hr"), list(k = i*12, p = (i+1)*12));
text(0, 0, range.name)
# with a variable, we can use for iteration to produce multiple tags
##############################

Plot in details

When we plot, we can assign the plot title, label.
We can also choose to show axis or not. On axis, we can set notations' letter size by cex.axis, direction by las.


Sample codes:

image.data = matrix(rnorm(60), nrow=10, ncol=6)
windows()
par(mar=c(5, 5, 4, 1))
image(x=1:nrow(image.data), y=1:ncol(image.data), image.data, axes = FALSE, xlab="X", ylab = "Y", main = "SAMPLE")
par(las=2) #control word's direction on plot, las = 1 or 2
namess1 = paste("row", 1:ncol(image.data), sep="-")
axis(2, 1:ncol(image.data), namess1, cex.axis=0.9, tick = FALSE)
par(las=2)
namess2 = paste("col", 1:nrow(image.data), sep="-")
axis(1, 1:nrow(image.data), namess2, cex.axis=0.9, tick = TRUE)
box()



windows()
par(mar=c(5, 5, 4, 1))
boxplot(data.frame(image.data), axes = FALSE, xlab="cols", ylab = "Distribution", main = "SAMPLE")
axis(2, tick = TRUE)
par(las=1)
namess2 = paste("col", 1:nrow(image.data), sep="-")
axis(1, 1:nrow(image.data), namess2, cex.axis=0.9, tick = TRUE)
box();