scatterplot3d, written by Uwe Ligges. Basic usage:
scatterplot3d(x,y,z)
scatterplot3d <-
function(x, y = NULL, z = NULL, color = par("col"), pch = NULL,
main = NULL, sub = NULL, xlim = NULL, ylim = NULL, zlim = NULL,
xlab = NULL, ylab = NULL, zlab = NULL, scale.y = 1, angle = 40,
axis = TRUE, tick.marks = TRUE, label.tick.marks = TRUE,
x.ticklabs = NULL, y.ticklabs = NULL, z.ticklabs = NULL,
y.margin.add = 0, grid = TRUE, box = TRUE, lab = par("lab"),
lab.z = mean(lab[1:2]), type = par("type"), highlight.3d = FALSE,
mar = c(5, 3, 4, 3) + 0.1, col.axis = par("col.axis"),
col.grid= "grey", col.lab= par("col.lab"), cex.symbols = par("cex"),
cex.axis = par("cex.axis"), cex.lab = 0.8 * par("cex.lab"),
font.axis = par("font.axis"), font.lab = par("font.lab"),
lty.axis = par("lty"), lty.grid = par("lty"), log = "", ...)
# log not yet implemented
{
## scatterplot3d, 0.3-11, 24.05.2002,
## Uwe Ligges ,
## http://www.statistik.uni-dortmund.de/leute/ligges.htm
##
## For MANY ideas and improvements thanks to Martin Maechler!!!
## Parts of the help files are stolen from the standard plotting functions in R.
mem.par <- par(mar = mar)
x.scal <- y.scal <- z.scal <- 1
xlabel <- if (!missing(x)) deparse(substitute(x))
ylabel <- if (!missing(y)) deparse(substitute(y))
zlabel <- if (!missing(z)) deparse(substitute(z))
## verification, init, ...
if(highlight.3d && !missing(color))
warning(message = "color is ignored when highlight.3d = TRUE")
if(length(x) < 2) stop("Minimal required length of x is 2 !")
## color as part of `x' (data.frame or list):
if(!is.null(d <- dim(x)) && (length(d) == 2) && (d[2] >= 4))
color <- x[,4]
else if(is.list(x) && !is.null(x$color))
color <- x$color
## convert 'anything' -> vector
xyz <- xyz.coords(x=x, y=y, z=z, xlab=xlabel, ylab=ylabel, zlab=zlabel,
log=log)
if(is.null(xlab)) { xlab <- xyz$xlab; if(is.null(xlab)) xlab <- "" }
if(is.null(ylab)) { ylab <- xyz$ylab; if(is.null(ylab)) ylab <- "" }
if(is.null(zlab)) { zlab <- xyz$zlab; if(is.null(zlab)) zlab <- "" }
if(length(color) == 1)
color <- rep(color, length(xyz$x))
else if(length(color) != length(xyz$x))
stop("length(color) must be equal length(x) or 1 !")
angle <- (angle %% 360) / 90
yz.f <- scale.y * abs(if(angle < 1) angle else if(angle > 3) angle - 4 else 2 - angle)
yx.f <- scale.y * (if(angle < 2) 1 - angle else angle - 3)
if(angle > 2) { ## switch y and x axis to ensure righthand oriented coord.
temp <- xyz$x; xyz$x <- xyz$y; xyz$y <- temp
temp <- xlab; xlab <- ylab; ylab <- temp
temp <- xlim; xlim <- ylim; ylim <- temp
}
angle.1 <- ifelse((2 > angle && angle > 1) || angle > 3, TRUE, FALSE)
angle.2 <- ifelse(1 > angle || angle > 3, FALSE, TRUE)
dat <- cbind(as.data.frame(xyz[c("x","y","z")]), col = color)
## xlim, ylim, zlim -- select the points inside the limits
if(!is.null(xlim)) {
xlim <- range(xlim)
dat <- dat[ xlim[1] <= dat$x & dat$x <= xlim[2] , , drop = FALSE]
}
if(!is.null(ylim)) {
ylim <- range(ylim)
dat <- dat[ ylim[1] <= dat$y & dat$y <= ylim[2] , , drop = FALSE]
}
if(!is.null(zlim)) {
zlim <- range(zlim)
dat <- dat[ zlim[1] <= dat$z & dat$z <= zlim[2] , , drop = FALSE]
}
n <- nrow(dat)
if(n < 1) stop("No data left within (x|y|z)lim")
y.range <- range(dat$y[is.finite(dat$y)])
if(all(diff(y.range) == 0))
stop("All points have the same Y-value! Use 2D-plot!")
### 3D-highlighting / colors / sort by y
if(type == "p" || type == "h") {
y.ord <- rev(order(dat$y))
dat <- dat[y.ord, ]
if(length(pch) > 1)
if(length(pch) != length(y.ord))
stop("length(pch) must be equal length(x) or 1 !")
else pch <- pch[y.ord]
if(highlight.3d)
dat$col <- rgb((1:n / n) * (y.range[2] - dat$y) / diff(y.range), g=0, b=0)
}
### optim. axis scaling
p.lab <- par("lab")
## Y
y.range <- range(dat$y, ylim)
y.prty <- pretty(y.range, n = lab[2],
min.n = max(1, min(.5 * lab[2], p.lab[2])))
y.scal <- round(diff(y.prty[1:2]), digits = 12)
y.add <- min(y.prty)
dat$y <- (dat$y - y.add) / y.scal
y.max <- (max(y.prty) - y.add) / y.scal
if(!is.null(ylim)) y.max <- max(y.max, ceiling((ylim[2] - y.add) / y.scal))
if(angle > 2) dat$y <- y.max - dat$y ## turn y-values around
## X
x.range <- range(dat$x[is.finite(dat$x)], xlim)
if(all(diff(x.range) == 0))
stop("All points have the same X-value! Use 2D-plot!")
x.prty <- pretty(x.range, n = lab[1],
min.n = max(1, min(.5 * lab[1], p.lab[1])))
x.scal <- round(diff(x.prty[1:2]), digits = 12)
dat$x <- dat$x / x.scal
x.range <- range(x.prty) / x.scal
x.max <- ceiling(x.range[2])
x.min <- floor(x.range[1])
if(!is.null(xlim)) {
x.max <- max(x.max, ceiling(xlim[2] / x.scal))
x.min <- min(x.min, floor(xlim[1] / x.scal))
}
x.range <- range(x.min, x.max)
## Z
z.range <- range(dat$z[is.finite(dat$z)], zlim)
if(all(diff(z.range) == 0))
stop("All points have the same Z-value! Use 2D-plot!")
z.prty <- pretty(z.range, n = lab.z,
min.n = max(1, min(.5 * lab.z, p.lab[2])))
z.scal <- round(diff(z.prty[1:2]), digits = 12)
dat$z <- dat$z / z.scal
z.range <- range(z.prty) / z.scal
z.max <- ceiling(z.range[2])
z.min <- floor(z.range[1])
if(!is.null(zlim)) {
z.max <- max(z.max, ceiling(zlim[2] / z.scal))
z.min <- min(z.min, floor(zlim[1] / z.scal))
}
z.range <- range(z.min, z.max)
### init graphics
plot.new()
if(angle.2) {x1 <- x.min + yx.f * y.max; x2 <- x.max}
else {x1 <- x.min; x2 <- x.max + yx.f * y.max}
plot.window(c(x1, x2), c(z.min, z.max + yz.f * y.max))
temp <- strwidth(as.character(y.scal * y.max + round(y.add, 0)),
cex = cex.lab/par("cex"))
if(angle.2) x1 <- x1 - temp - y.margin.add
else x2 <- x2 + temp + y.margin.add
plot.window(c(x1, x2), c(z.min, z.max + yz.f * y.max))
if(angle > 2) par("usr" = par("usr")[c(2, 1, 3:4)])
title(main, sub, ...)
### draw axis, tick marks, labels, grid, ...
if(grid) {
## X
i <- x.min:x.max
segments(i, z.min, i + (yx.f * y.max), yz.f * y.max + z.min,
col = col.grid, lty = lty.grid)
## Y
i <- 0:y.max
segments(x.min + (i * yx.f), i * yz.f + z.min,
x.max + (i * yx.f), i * yz.f + z.min,
col = col.grid, lty = lty.grid)
}
if(tick.marks && axis) { ## tick marks
xtl <- (z.max - z.min) * (tcl <- -par("tcl")) / 50
ztl <- (x.max - x.min) * tcl / 50
## Y
i <- 0:y.max
temp <- ifelse(angle.2, x.min, x.max)
segments(yx.f * i - ztl + temp, yz.f * i + z.min,
yx.f * i + ztl + temp, yz.f * i + z.min,
col=col.axis, lty=lty.axis)
## X
i <- x.min:x.max
segments(i, -xtl + z.min, i, xtl + z.min, col=col.axis, lty=lty.axis)
## Z
i <- z.min:z.max
temp <- ifelse(angle.2, x.max, x.min)
segments(-ztl + temp, i, ztl + temp, i, col=col.axis, lty=lty.axis)
if(label.tick.marks) { ## label tick marks
las <- par("las")
mytext <- function(labels, side, at, ...)
mtext(text = labels, side = side, at = at, line = -.5,
col=col.lab, cex=cex.lab, font=font.lab, ...)
## X
j <- subset(temp <- pretty(x.range, n = lab[1]),
temp <= x.range[2] & temp >= x.range[1])
if(is.null(x.ticklabs))
x.ticklabs <- j * x.scal
mytext(x.ticklabs, side = 1, at = j)
## Z
j <- subset(temp <- pretty(z.range, n = lab.z),
temp <= z.range[2] & temp >= z.range[1])
if(is.null(z.ticklabs))
z.ticklabs <- j * z.scal
mytext(z.ticklabs, side = ifelse(angle.1, 4, 2), at = j,
adj = ifelse((0 < las) && (las < 3), 1, NA))
## Y
j <- subset(temp <- pretty(c(0, y.max), n = lab[2]),
temp <= y.max & temp >= 0)
temp <- if(angle > 2) rev(j) else j ## turn y-labels around
if(is.null(y.ticklabs))
y.ticklabs <- y.scal * temp + round(y.add, 0)
else if (angle > 2)
y.ticklabs <- rev(y.ticklabs)
text(j * yx.f + ifelse(angle.2, x.min, x.max),
j * yz.f + z.min, y.ticklabs,
pos = ifelse(angle.1, 2, 4), offset = 1,
col=col.lab, cex = cex.lab/par("cex"), font=font.lab)
}
}
if(axis) { ## axis and labels
mytext <-
function(lab, side = side, at = at, ...)
mtext(lab, side = side, at = at, col = col.lab,
cex = cex.axis, font = font.axis, las = 0, ...)
## X
lines(c(x.min, x.max), c(z.min, z.min), col = col.axis, lty = lty.axis)
mytext(xlab, line = 1.5, side = 1, at = mean(x.range))
## Y
lines(c(x.max, x.max + y.max * yx.f), c(z.min, y.max * yz.f + z.min),
col = col.axis, lty = lty.axis)
mytext(ylab, side = ifelse(angle.1, 2, 4), at = z.min + y.max * yz.f, line = .5)
## Z
lines(c(x.min, x.min), c(z.min, z.max), col = col.axis, lty = lty.axis)
mytext(zlab, line = 1.5, side = ifelse(angle.1, 4, 2), at = mean(z.range))
if(box) {
## X
temp <- yx.f * y.max
temp1 <- yz.f * y.max
lines(c(x.min + temp, x.max + temp),
c(z.min + temp1, z.min + temp1), col = col.axis, lty = lty.axis)
lines(c(x.min + temp, x.max + temp), c(temp1 + z.max, temp1 + z.max),
col = col.axis, lty = lty.axis)
## Y
temp <- c(0, y.max * yx.f)
temp1 <- c(0, y.max * yz.f)
lines(temp + x.min, temp1 + z.min, col = col.axis, lty = lty.axis)
lines(temp + x.min, temp1 + z.max, col = col.axis, lty = lty.axis)
## Z
temp <- yx.f * y.max
temp1 <- yz.f * y.max
lines(c(temp + x.min, temp + x.min), c(z.min + temp1, z.max + temp1),
col = col.axis, lty = lty.axis)
lines(c(x.max + temp, x.max + temp), c(z.min + temp1, z.max + temp1),
col = col.axis, lty = lty.axis)
}
}
### plot points
x <- dat$x + (dat$y * yx.f)
z <- dat$z + (dat$y * yz.f)
col <- as.character(dat$col)
if(type == "h") {
z2 <- dat$y * yz.f + z.min
segments(x, z, x, z2, col = col, cex = cex.symbols, ...)
points(x, z, type = "p", col = col, pch = pch, cex = cex.symbols, ...)
}
else points(x, z, type = type, col = col, pch = pch, cex = cex.symbols, ...)
### box-lines in front of points (overlay)
if(axis && box) {
lines(c(x.min, x.max), c(z.max, z.max),
col = col.axis, lty = lty.axis)
lines(c(0, y.max * yx.f) + x.max, c(0, y.max * yz.f) + z.max,
col = col.axis, lty = lty.axis)
lines(c(x.max, x.max), c(z.min, z.max),
col = col.axis, lty = lty.axis)
}
par(mem.par)
### Return Function Object
ob <- ls() ## remove all unused objects from the result's enviroment:
rm(list = ob[!ob %in% c("x.scal", "y.scal", "z.scal", "yx.f",
"yz.f", "y.add", "z.min", "z.max", "x.min", "x.max", "y.max")])
rm(ob)
invisible(list(
xyz.convert = function(x, y=NULL, z=NULL) {
xyz <- xyz.coords(x, y, z)
y <- (xyz$y - y.add) / y.scal
return(x = xyz$x / x.scal + yx.f * y,
y = xyz$z / z.scal + yz.f * y)
},
points3d = function(x, y = NULL, z = NULL, type = "p", ...) {
xyz <- xyz.coords(x, y, z)
y2 <- (xyz$y - y.add) / y.scal
x <- xyz$x / x.scal + yx.f * y2
y <- xyz$z / z.scal + yz.f * y2
if(type == "h") {
y2 <- z.min + yz.f * y2
segments(x, y, x, y2, ...)
points(x, y, type = "p", ...)
}
else points(x, y, type = type, ...)
},
plane3d = function(Intercept, x.coef = NULL, y.coef = NULL, lty = "dashed", ...){
if(!is.null(coef(Intercept))) Intercept <- coef(Intercept)
if(is.null(x.coef) && length(Intercept) == 3){
x.coef <- Intercept[2]
y.coef <- Intercept[3]
Intercept <- Intercept[1]
}
x <- x.min:x.max
x.coef <- x.coef * x.scal
z1 <- (Intercept + x * x.coef + y.add * y.coef) / z.scal
z2 <- (Intercept + x * x.coef +
(y.max * y.scal + y.add) * y.coef) / z.scal
segments(x, z1, x + y.max * yx.f, z2 + yz.f * y.max, lty = lty, ...)
y <- 0:y.max
y.coef <- (y * y.scal + y.add) * y.coef
z1 <- (Intercept + x.min * x.coef + y.coef) / z.scal
z2 <- (Intercept + x.max * x.coef + y.coef) / z.scal
segments(x.min + y * yx.f, z1 + y * yz.f,
x.max + y * yx.f, z2 + y * yz.f, lty = lty, ...)
},
box3d = function(...){
lines(c(x.min, x.max), c(z.max, z.max), ...)
lines(c(0, y.max * yx.f) + x.max, c(0, y.max * yz.f) + z.max, ...)
lines(c(x.max, x.max), c(z.min, z.max), ...)
lines(c(x.min, x.max), c(z.min, z.min), ...)
}
))
}
z <- seq(-10,10,.01)
x <- cos(z)
y <- sin(z)
scatterplot3d(x,y,z, highlight.3d=TRUE, col.axis="blue",
col.grid="lightblue",pch=20)
Income <- state.x77[,2]
Illit <- state.x77[,3]
LifeExp <- state.x77[,4]
scatterplot3d(Income, Illit, LifeExp)
scatterplot3d(Income, Illit, LifeExp,type="h")
scatterplot3d(Income, Illit, LifeExp,
highlight.3d=TRUE,type="h",pch=15)