You are here

Three dimensional plots using rgl package

We all know that R can do amazing things including 3 dimensional plots. scatterplot3d is probably the most popular package for doing this. But a few days ago I got introduced with rgl package which can do 3 dimensional graphs with some added advantages like we can rotate the plot using mouse, zoom in or out using the mouse scroll wheel and even can play beautiful animations.
The data I have used for all the examples here is-

    X    Y       Z
 2276 1902   156.6
 2925 2138  -274.4
 2209 1777    57.7
 2712 1603   269.4
 3128 1545   391.4
 .... ....   .....
13526 8575  2178.8
12268 8063  3291.4
12006 7747  4402.9
10763 8752  2533.3
11927 9932   156.6
 

Here, X,Y and Z are monthly import,export and first differenced net foreign asset of Bangladesh from July 1998 to July 2010 respectively. That means there are a total of 133 observations in each vector. Here is the command for doing a scatter plot in rgl package-

plot3d(x=X, y=Y, z=Z, type="p", col="red", xlab="X", ylab="Y", zlab="Z", 
size=5, lwd=15, box=F)

And the plot is (we can rotate the plot as our will, here just two different views are given)-
mei112.png
mei113.png
We can save the screen shot of the plot as png file by the command-

rgl.snapshot(filename="mei11.png",fmt="png")

The plot can be saved as other format like pdf and eps by rgl.postscript command.
In either cases the saved file will no longer be rotatable.
Now, we can make the plot a bit more beautiful by putting type="s" in plot3d function. That means the points will be represented by spheres. The code is-

plot3d(x=X, y=Y, z=Z, type="s", col="red", xlab="X", ylab="Y", zlab="Z",
size=5, radius=200, box=F)

mei114.png

If we leave the radius argument as default and put col=rainbow(n) the plot becomes very interesting to see. The command will be-

plot3d(x=X, y=Y, z=Z, type="s", col=rainbow(n=10), xlab="X", ylab="Y",
zlab="Z", size=5, box=F)

And the output of this command is the main figure of this article, have a look again.
This plot gives us no added advantage, in fact it is very difficult to understand, but I can not resist myself (because the plot is so beautiful now!).

Now, let’s do some fancy animations. Let’s rotate the plot automatically so that we can relax and see every aspect of our three dimensional plot.
Here are the commands-

M <- par3d("userMatrix")
play3d( par3dinterp( userMatrix=list(M,rotate3d(M, angle=pi,x=1 , y=0,z= 0) ) ), duration=10)

This command will rotate the plot for 10 seconds in the X axis at 180 degrees.
If you want to rotate it in Y axis, just put x=0 and y=1. And by changing the values of angle, x, y and z we can make the plot rotate in every possible way.

Now let’s do something more interesting- let’s write a simple function that would make the plot rotate automatically in several ways for a given number of times, so that we can have very well and details look of it.

The function is-

plot3d(x=X,y=Y,z=Z,type="s",col=rainbow(n=10),xlab="X",ylab="Y", zlab="Z",size=5,lwd=15,box=F)
M <- par3d("userMatrix")
st<-function(n,t)  {
            for (i in 1:n) {
              play3d( par3dinterp( userMatrix=list(M,rotate3d(M, pi, 1, 0, 0) ) ), duration=t )
             play3d( par3dinterp( userMatrix=list(M,rotate3d(M, pi, 0, 1, 1) ) ), duration=t)
             play3d( par3dinterp( userMatrix=list(M,rotate3d(M, pi, 0, 1, 0) ) ), duration=t)
             play3d( par3dinterp( userMatrix=list(M,rotate3d(M, pi, 1, 0, 1) ) ), duration=t)
             }
}

Now write-

st(20,2)

Watch and enjoy.

There are so many other interesting things in rgl package. Those who are interested in statistical programming (especially in R) please try it and share your knowledge and experience and help me to enhance mine.

Category: 
Tag: 

Comments

Good topic Tanvir. Thanks for posting.

If you add some reproducible data that would benefit average users like me. For example, in you code what is X, Y, Z? I assume that they are vectors, right?

In short, an article is most useful if it is self contained. An example with a real example helps understanding things better.

Another thing is, if you have a main figure for the article, don't show it again in the body of the article. This will save some bandwidth for people in Bangladesh. Instead, you can just refer to the figure above and give the code near the bottom.

Keep up the good work.

Thank you sir. I have tried to fix the pitfalls. The data has 133 rows, so I put just a few of them, is it OK?
Good joke sir ("....average users like me...."). It sounds funny when a very much capable and experienced R user (and statistician) like you call himself "average user".
Again, thank you sir for your very much encouraging and suggestive words.

Here is one particular example of 3d plots using not-so-fancy packages (basic and lattice):

A brief description of the problem:

gk-distribution is a broader class of distribution with parameters G (skewness) and K (kurtosis). In this family of distribution, G = 0 and K = 0 gives a special case - which is popularly known as normal distribution (skewness and kurtosis values are zero). Thus, to check the performance of a normality test such asPearson Chi-square test, a simple strategy would be to generate data from gk-distribution with various values of G and K, and then see whether according to the test, those data follows normality or not, in the entire parameter grid. Similar strategy would be also applicable to check the effect of sample size. As we know from theory, the greater the sample size, the better the power curves should be. Therefore, in this analysis, I generate data from gk-distribution using sample sizes 5, 10, 15 and 30. Also to produce the grid, I used values from -5 to +5 for G, and -0.5 to 5 for K. Since I am basically interested what happens near (0, 0) location, the grid is tighter there. Since this requires 3D plots to show the powers of the parameter grid, I used both basic and lattice plots. Eventually it can be shown that the more the sample size, the better the powers, since the rejection rate near (G, K) = (0, 0) is very low for higher sample sizes, and very high in all other locations. However, when the sample size gets smaller, the performance of the powers are not very encouraging in favor of this test.

Click on the hyperlinked code to see the output:

# load required packages
require(nortest)
require(lattice)

# Setting the sample size and number of simulations.
# If the "sampleSize.fix" value is changed and sourced from R,
# this program will automatically generate graphs
# for specifies sample size
sampleSize.fix = 30 # The other values of sample size are 5, 10, 15
simNum.fix = 100

# Function to generate gk distribution data
gk <- function(G, K, simNum, sampleSize){
rejectionCount <- rep(0,simNum)
for(i in 1:simNum){
z <- rnorm(sampleSize, 0, 1)
gkGenerated <- 0+1*z*(1+.83*(1-exp(-G*z))/(1+exp(-G*z)))*(1+z^2)**K
gkTest <- pearson.test(gkGenerated)
pearsonPvalue <- gkTest\$p.value
if(pearsonPvalue <= 0.05) {
rejectionCount[i] <- 1
}
else {
rejectionCount[i] <- 0
}
}
power <- mean(rejectionCount)
return(power)
}

# Making the grid with G (skewness) and K (kurtosis) parameters
G <- c(-5, -4, -3, -2, -1.5, -1, -0.75, -0.5, -0.25,
0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5)
K <- c(-0.5, -0.25, 0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5)
lg <- length(G)
lk <- length(K)
parameterGrid <- expand.grid(Ggrid = G, Kgrid = K)
parameterGrid <- as.data.frame(parameterGrid)
ggrid <- parameterGrid[[1]]
kgrid <- parameterGrid[[2]]

# Linking gk function with the parameter grid
gkPower <- function(x, simNum, sampleSize){
gk(ggrid[x], kgrid[x], simNum, sampleSize)
}

gkPowerStore <- double(length = (lg*lk))
# storing the powers of analysis
for (i in 1:(lg*lk)) {
gkPowerStore[i] <- gkPower(i, simNum = simNum.fix,
sampleSize = sampleSize.fix)
cat("iteration", i, "\n")
}
# converting the vector into matrix for convenience of making the plots
powerMatrix <- matrix(gkPowerStore, nrow=lg, ncol=lk, byrow=F)

# producing basic plots
x11()
persp(G, K, powerMatrix, theta = 30, phi = 30,
expand = 0.5, d = 1, ltheta = 120,
shade = 0.5, ticktype = "detailed",
xlab = "Parameter G", ylab = "Parameter K",
zlab = "Power", col = "lightblue",
main = "Basic Graph", sub = paste("Sample size =", sampleSize.fix))
dev.copy2pdf(file = paste("basic", sampleSize.fix, ".pdf", sep = ""),
onefile = TRUE, paper = "a4",family = "Helvetica",
pointsize=1, fonts = "Helvetica")
dev.copy2eps(file = paste("basic", sampleSize.fix, ".eps", sep = ""),
onefile = TRUE, paper = "a4",family = "Helvetica",
pointsize=1, print.it = FALSE, fonts ="Helvetica")
dev.off()

# producing lattice plots
x11()
wireframe(powerMatrix~parameterGrid\$Kgrid*parameterGrid\$Ggrid,
main = "Lattice Graph", xlab = "Parameter K", ylab = "Parameter G",
zlab = "Power", shade = T,
sub = paste("Sample size =", sampleSize.fix))
dev.copy2pdf(file = paste("lattice", sampleSize.fix, ".pdf", sep = ""),
onefile = TRUE, paper = "a4",family = "Helvetica",
pointsize=1, fonts = "Helvetica")
dev.copy2eps(file = paste("lattice", sampleSize.fix, ".eps", sep = ""),
onefile = TRUE, paper = "a4",family = "Helvetica", pointsize=1,
print.it = FALSE, fonts ="Helvetica")
dev.off()

Warning: Running this piece of code will generate 4 graphs in your R working directory (see getwd).

Thanks Ehsan vaia, you have given something great. I am so honored and glad that you have written my post and also commented on it. But the function will really take sometime for me to digest (I am a very much beginner).
I have given the codes for rotating the plot automatically, your valuable suggestion,comment and more such functions are desired....

Thank you again vaia.

Don't read too much into it, ignore the data generating procedure (i just gave something extra, if someone wants to get the whole perspective of the idea): just see the format of powerMatrix (print it in the console), and then see the use of persp (basic) and wireframe (lattice) :)

Thanks again vaia, now I have started to gain some courage to try the function.... In my last comment I made a mistake; I wrote "you have written..." it would be "you have read my post and also commented on it" instead.
Hope we all will be able to learn many other amazing things from you.
Good wishes vaia.

I just upgraded to R version 2.12.1 and rgl version 0.92.798. Unlike the previous versions, rgl.snapshot() no longer works. In other words, except for a screen save, there's no way that I can find to save the image/movie.

Do you know of a work-around?

Is there a way to inform those that are maintaining the package?

George

It is unfortunate that the most recent version of rgl package does not support rgl.snapshot command.But you can use rgl.postscript command to save the 3d plot in ps,eps or pdf format. The command is-

rgl.postscript(filename="mei11.pdf",fmt="pdf")

Use gimp or photoshop(in windows) to convert the pdf(or ps/eps) plot to your favourite format(like png).

Hope this will help you.

There are now several of us in the same boat. I just got another panic e-mail.

It appears to be rgl version that is the problem. And, the postscript route may or may not work. I have several runs where it displays the right image, but the saved .eps file gives the wrong shading and colors.

I haven't tested it yet, but one guy is getting huge files. This is a big problem.

Yes the saved eps file is not of good quality. Probably it would be better if you first save the plot in pdf format and then convert it into png(or any other format). I am not happy with the saved pdf image also, but it is certainly better than the eps. Yes, this is a big problem. Hope the next version will support snapshot again.

Good day!
I've been surfing the net for long hours to look for R command in plotting 2D and 3Dimensional solution space... Fortunately, i found this page..
I just want to ask if there is R command to plot the data points and make use of symbols as group indicator.. to come up with a conclusion that a points of the same symbol belongs to the same group.. something like that.. it find me hard to do so..
Any suggestions please.. any help will do and very much appreciated.. THANKS!

Let, we want to plot 2 vectors,x and y, on a two dimensional space,and we want different symbols and colors for each of the situations-
(1) x is less than y
(2) x=y
(3) x>y

Let the vectors are-
x<-c(1,2,1,5,4,6,7,3,4,2,1)
y<-c(0,2,3,5,7,2,0,1,2,3,4)

Now, try the codes-

p<-function(x,y) {
   p<-c()
   n<-length(x)
   for (i in 1:n) {
   if (x[i]<y[i]) p[i]<-1
   else if (x[i]==y[i]) p[i]<-2
   else p[i]<-3
   }
   return(p)
}

And finally-
plot(x,y,pch=p(x,y),col=p(x,y))

Is this the thing you are trying to do?

Add new comment