ContextBase Logo



This is a demonstration of 3d plotting in the R programming language



Required Packages

setwd("C:/Users/Administrator/Dropbox/Programming/DataVisualization")

# install.packages("plot3D")
library(plot3D)

# install.packages("rgl")
library(rgl)

# https://cran.r-project.org/web/packages/scatterplot3d/vignettes/s3d.pdf
# install.packages("scatterplot3d")
library(scatterplot3d)



Basic 3D plotting in R

fun <- function(x, y) { (x+y)/2 }
x <- y <- seq(1, 3, length=20)
z <- outer(x,y,fun)
persp(x, y, z)

x <- seq(0, 20, length=30)
y <- x
fun <- function(x, y) { (x+y)/2 }
z <- outer(x, y, fun)
persp(x, y, z, theta=50, phi=20, expand=0.5, col="blue")

x <- seq(-10, 10, length=30)
y <- x
persp(x, y, z, theta=50, phi=20, expand=0.5, col="green",
      ltheta=60, shade=0.75, ticktype="detailed",
      xlab="X", ylab ="Y", zlab=" r ")

y <- x <- seq(3,6,length=30)
fun <- function(x,y) { return(dnorm(x)/dnorm(y)) }
z <- outer(x,y,fun)
persp(x,y,z)



3D Plot of Half of a Torus

if (!require("plot3D")) {install.packages("plot3D"); require("plot3D")}
library(plot3D)
par(mar = c(2, 2, 2, 2))
par(mfrow = c(1, 1))
R <- 3
r <- 2
x <- seq(0, 2*pi,length.out=50)
y <- seq(0, pi,length.out=50)
M <- mesh(x, y)
 
alpha <- M$x
beta <- M$y
 
surf3D(x = (R + r*cos(alpha)) * cos(beta),
       y = (R + r*cos(alpha)) * sin(beta),
       z = r * sin(alpha), colkey=FALSE, bty="b2",
       main="Half of a Torus")



Scatterplot3D

# https://cran.r-project.org/web/packages/scatterplot3d/vignettes/s3d.pdf
if (!require("scatterplot3d")) {install.packages("scatterplot3d"); require("scatterplot3d")}
library(scatterplot3d)
z <- seq(-10, 10, 0.01)
x <- cos(z)
y <- sin(z)
scatterplot3d(x, y, z, highlight.3d = TRUE, col.axis = "blue",
col.grid = "lightblue", main = "Helix", pch = 20)

data(trees)
s3d <- scatterplot3d(trees, type = "h", color = "blue",
angle = 55, scale.y = 0.7, pch = 16, main = "Adding elements")
my.lm <- lm(trees$Volume ~ trees$Girth + trees$Height)
s3d$plane3d(my.lm)
s3d$points3d(seq(10, 20, 2), seq(85, 60, -5), seq(60, 10, -10),
col = "red", type = "h", pch = 8)

cubedraw <- function(res3d, min = 0, max = 255, cex = 2)
{
cube01 <- rbind(0,c(1,0,0),c(1,1,0),1,c(0,1,1),c(0,0,1),c(1,0,1),
c(1,0,0),c(1,0,1),1,c(1,1,0),
c(0,1,0),c(0,1,1), c(0,1,0),0)
cub <- min + (max-min)* cube01
res3d$points3d(cub[ 1:11,], cex = cex, type = 'b', lty = 1)
res3d$points3d(cub[11:15,], cex = cex, type = 'b', lty = 3)
}
crgb <- t(col2rgb(cc <- colors()))
rr <- scatterplot3d(crgb, color = cc, box = FALSE, angle = 24)
cubedraw(rr)

Rrb <- t(col2rgb(rbc <- rainbow(201)))
rR <- scatterplot3d(Rrb, color = rbc, box = FALSE, angle = 24)
cubedraw(rR)
rR$points3d(Rrb, col = rbc, pch = 16)



Population Data 3D Plot

data(USPersonalExpenditure)
data(uspop)

expendData <- USPersonalExpenditure[1:4,]
expendData <- cbind(expendData[,1], expendData[,3], expendData[,5])
colnames(expendData) <- c("1940", "1950", "1960")
populationData <- uspop
populationData <- data.frame(populationData)
expendData <- rbind(expendData, populationData[16:18,])
rownames(expendData) <- c("Food and Tobacco", "Household Operation",
                          "Medical and Health", "Personal Care",
                          "Population")

library(knitr)
kable(expendData, caption="USA Personal Expenditure Data")
USA Personal Expenditure Data
1940 1950 1960
Food and Tobacco 22.20 59.60 86.8
Household Operation 10.50 29.00 46.2
Medical and Health 3.53 9.71 21.1
Personal Care 1.04 2.45 5.4
Population 131.70 151.30 179.3 
plot(USPersonalExpenditure)

plot(uspop)

persp(USPersonalExpenditure, theta=40, phi=10, expand=0.5,
      col="green", ltheta=60, shade=0.75, ticktype="detailed",
      xlab="", ylab="", zlab="")



Perspective Box

perspbox(z = volcano, bty = "g", 
         d = 2, main  = "3D Coordinates")



Paraboloid

# inversion function
# f1 <- function(x,y){
#   x^2 + y^2
# }

f1 <- function(x,y){
  -x^2 - y^2 + 10
}
x <- seq(-150, 150,length = 100)
y <- x
z <- outer(x,y,f1)

persp(x, y, z, theta=-10, phi=20, expand=1,
      col = "blue", shade = 0.75, border = NA,
      ticktype = "detailed", main = "Pltfiv")



Surf3D Paraboloid

x <- seq(-8, 8 , length.out=90)
y <- seq(-8, 8, length.out=90)
M <- mesh(x, y)

alpha <- M$x
beta <- M$y

surf3D(x = M$x, y = M$y, z = -(alpha-1)^2 - beta^2+20,
       theta=-10, phi=20, colvar=-(alpha-1)^2 - beta^2+20,
       colkey=TRUE, bty="b2", main ="Paraboloid")



rgl Paraboloid With Intersecting Plane

x <- seq(-8, 8 , length.out=90)
y <- seq(-8, 8, length.out=90)
M <- mesh(x, y)

alpha <- M$x
beta <- M$y

plot3d(x = M$x , y = M$y, z = -(alpha-1)^2 - beta^2+20,
       theta=-10, phi=20,
       main ="Paraboloid with Intersecting Plane")
planes3d(4, -4, 4, 0, col = 'red', alpha = 0.6)

rgl Paraboloid



Paraboloid with Intersecting Plane

library(rgl)
x = seq(-4, 4, 0.2)
y = seq(-4, 4, 0.2)
z = matrix(data=NA, nrow=length(x), ncol=length(x))
# fill the z matrix
sigma = 1.0; mux = 0.0; muy = 0.0; A = 1.0

for(i in 1:length(x)) {
  for(j in 1:length(y))
  {z[i,j] =  A * (1/(2*pi*sigma^2)) *
    exp( -((x[i]-mux)^2 + (y[j]-muy)^2)/(2*sigma^2) ) 
  }}

persp3d(x, y, z, nticks=5, ticktype="detailed",
        zlim = c(0,0.2), col="cyan", xlab="X",
        ylab="Y", zlab="Z", space = 0.1, axes = TRUE,
        main ="Paraboloid with Intersecting Plane")

# x and y = corners of the intersecting plane, z = height
sqdf <- data.frame(x=c(1,1,1,1,1),
                   y=c(-3,3,3,-3,-3),
                   z=c(0.125,0.125,0,0,0.125))

# draw the intersecting plane, 
#    the "add=T" parameter appends it to the previous 3d-plot
#    the coord paramter tells it what two planes to use when 
#    tesselating the polygon into triangles

polygon3d(sqdf$x, sqdf$y, sqdf$z, coord=2:3, alpha=0.5,
          color="purple", add=TRUE)
view3d(theta = 20, phi = -60)

Paraboloid with Intersecting Plane



Partial Derivation

z <- volcano[40:87, 30:50]
x <- 10 * (1:nrow(z))
y <- 10 * (1:ncol(z))
z0 <- min(z) - 20
z <- rbind(z0, cbind(z0, z, z0), z0)
x <- c(min(x) - 1e-10, x, max(x) + 1e-10)
y <- c(min(y) - 1e-10, y, max(y) + 1e-10)

fill <- matrix("green3", nrow = nrow(z)-1, ncol = ncol(z)-1)
fill[ , i2 <- c(1,ncol(fill))] <- "gray"
fill[i1 <- c(1,nrow(fill)) , ] <- "gray"
par(bg = "lightblue")

persp3d(x, y, z, col = fill,
        main = "Partial Derivation")

# x and y = corners of the intersecting plane, z = height
sqdf <- data.frame(x=c(250,250,250,250,250),
                   y=c(0,200,200,0,0),
                   z=c(150,150,80,80,150))

polygon3d(sqdf$x, sqdf$y, sqdf$z, coord=2:3, alpha=0.5,
          color="purple", add=TRUE)
view3d(theta = 20, phi = -60)

Partial Derivation



3D Plot Animation

library(animation)

saveGIF({
  for(i in 1:150){
    x <- seq(-6 + (i * 0.05), 6 + (i * 0.05), length= 100)
    y <- x
    f <- function(x, y) { sin(x) + cos(y) }
    z <- outer(x, y, f)
    persp(x, y, z, theta = 45 + (i * 0.5), phi = 35, expand = 0.4, col = "lightblue")
  }
}, interval = 0.1, ani.width = 550, ani.height = 550)

3D Plot Animation