## Other 3D Viewing Options

Other options for viewing a 3D scene include:

A stereogram

• presents each eye with an image from a slightly different viewing angle;
• the brain fuses the images into a 3D view in a process called stereopsis.

This approach is or was used

You can create your own viewer with paper or cardboard tubes, for example.

A stereo image for the soil data surface:

tp <- trellis.par.get()
trellis.par.set(theme=col.whitebg())
ocol <- trellis.par.get("axis.line")\$col
oclip <- trellis.par.get("clip")\$panel
trellis.par.set(list(axis.line = list(col = "transparent"),
clip = list(panel = FALSE)))
print(wireframe(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9,
aspect=diff(range(soi.grid\$n))/diff(range(soi.grid\$e)),
xlab = "Easting (km)",
screen = list(z = 40, x = -70, y = 3)),
split = c(1,1,2,1), more = TRUE)

print(wireframe(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9,
aspect=diff(range(soi.grid\$n))/diff(range(soi.grid\$e)),
xlab = "Easting (km)",
screen = list(z = 40, x = -70, y = 0)),
split = c(2,1,2,1))

trellis.par.set(list(axis.line = list(col = ocol), clip = list(panel = oclip)))
trellis.par.set(tp)

## Some Calibration Experiments

### Simple Linear Stucture

Linear structure with normally distributed data:

n <- 1000
d1 <- data.frame(x = rnorm(n), y = rnorm(n))
d1 <- mutate(d1, z = x + y)
splom(d1)

xyplot(z ~ x | equal.count(y), data = d1)

xyplot(z ~ x | equal.count(y, 12), data = d1)

xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d1)

splom(d1, group = abs(x) < 0.1)

clear3d()
points3d(d1)

Linear structure with uniformly distributed data:

d2 <- data.frame(x = runif(n, -1, 1), y = runif(n, -1, 1))
d2 <- mutate(d2, z = x + y)
splom(d2)

xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d2)

splom(d2, group = abs(x) < 0.1)

clear3d()
points3d(d2)

### Nonlinearity and Noise

Normal data with non-linearity and noise:

d3 <- mutate(d1, z = cos(x) + cos(y) + 0.5 * (x + y) + rnorm(x, 0, .2))
splom(d3)

xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d3)

splom(d3, group = abs(x) < 0.1)

splom(d3, group = abs(x + 1) < 0.1)

Uniform data with non-linearity and noise:

d4 <- mutate(d2,  z = cos(x) + cos(y) + 0.5 * (x + y) + rnorm(x, 0, .1))
xyplot(z ~ x | equal.count(y, 12, overlap = -0.8), data = d4)

splom(d4, group = abs(x) < 0.1)

Full data in the background for context:

xyplot(z ~ x | equal.count(y, 12, overlap = -0.8), data = d4,
panel = function(x, y, ...) {
panel.xyplot(d4\$x, d4\$z, col = "grey")
panel.xyplot(x, y,...)
})

clear3d()
points3d(d4)

### More Variables

Two-dimensional views, conditioning on two variables:

n <- 10000
d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
d <- mutate(d,
x4 = as.numeric(scale(x1 + cos(x2) + sin(x3) + rnorm(n, 0, 0.1))))
splom(d)

splom(~d, data = d, group = abs(x1) < 0.1)

splom(~d, data = d, group = abs(x1) < 0.1 & abs(x2) < 0.1)

xyplot(x4 ~ x3 |
equal.count(x1, 4, overlap = -0.8) +
equal.count(x2, 4, overlap = -0.8),
data = d, asp = 1)

Three-dimensional views, conditioning on one variable:

dd <- filter(d, abs(x1) < 0.1)[2:4]
clear3d()
points3d(dd)

clear3d()

clear3d()
spheres3d(dd, radius = 0.05, col = ifelse(abs(dd\$x2) < 0.1, "red", "blue"))

dd1 <- filter(d, abs(x1 - 1) < 0.1)[2:4]
dd2 <- filter(d, abs(x1 + 1) < 0.1)[2:4]
clear3d()
spheres3d(dd1, radius = 0.05, color = rep("red", ncol(dd1)))
spheres3d(dd2, radius = 0.05, color = rep("blue", ncol(dd2)))

ddx <- mutate(d,
v = ifelse(abs(x1 - 1) < 0.1,
1,
ifelse(abs(x1 + 1) < 0.1, 2, NA)))
ddx <- filter(ddx, ! is.na(v))
cloud(x4 ~ x2 * x3, group = v, data = ddx, auto.key = TRUE)

cloud(x4 ~ x2 * x3 | v, data = ddx)

## Interactive Explorations

### Tk Sliders Using tkrplot

A simple animation function using the tkrplot` package:

library(tkrplot)

tkranim <- function(fun, v, lo, hi, res = 0.1, title, label) {
tt <- tktoplevel()
if (! missing(title))
tkwm.title(tt, title)
draw <- function() fun(v)
img <-tkrplot(tt, draw)
tv <- tclVar(init = v)
f<-function(...) {
vv <- as.numeric(tclvalue(tv))
if (vv != v) {
v <<- vv
tkrreplot(img)
}
}
s <- tkscale(tt, command=f, from=lo, to=hi, variable=tv,
showvalue=FALSE, resolution=res, orient="horiz")
if (! missing(label))
tkpack(img, tklabel(tt, text=label), s)
else tkpack(img,s)
}

Interactive conditioning for one of the calibration examples:

fd4 <- function(v)
print(splom(~ d4, group = abs(x - v) < 0.1, data = d4))

tkranim(fd4, 0, -1, 1, 0.05)

Interactive slicing for the soil resistivity surface:

f <- function(i) {
e <- eastseq[i]
idx <- which(sf\$easting == e)
r <- range(sf\$fit)
print(xyplot(fit ~ northing, data = sf[idx,],
ylim = r, type = "l", main = sprintf("Easting = %f", e)))
}

tkranim(f, 1, 1, 85, 1)

Interactive conditioning for the soil resistivity data:

binwidth <- 0.1 ## bad idea to use a global variable
g <- function(i) {
e <- eastseq[i]
idx <- which(abs(soil\$easting - e) <= binwidth)
ry <- range(soil\$resistivity)
rx <- range(soil\$northing)
print(xyplot(resistivity ~ northing, data = soil[idx,],
xlim = rx, ylim = ry, main = sprintf("Easting = %f", e)))
}

tkranim(g, 1, 1, 85, 1)

### Shiny Sliders

A version of the animation function using the shiny package:

shinyAnim <- function(fun, v, lo, hi, step = NULL, label = "Frame") {
require(shiny)
shinyApp(
ui = fluidPage(
sidebarLayout(
sidebarPanel(sliderInput("n", label, lo, hi, v, step,
animate = animationOptions(100))),
mainPanel(plotOutput("plot"))
)
),
server = function(input, output, session) {
session\$onSessionEnded(function() shiny::stopApp())
output\$plot <- renderPlot(fun(input\$n))
}
)
}
shinyAnim(fd4, 0, -1, 1, 0.05)
shinyAnim(f, 1, 1, 85)
shinyAnim(g, 1, 1, 85, 1)

### Animations as Movies

To show how a plot changes with a single parameter it is possible to pre-compute the images and assemble them into a movie.

This is easy to do in Rmarkdown, but requires the ffmpeg software to be installed.

Currently, ffmpeg is not available on the CLAS Linux systems.

Linked plots and brushing were pioneered at Bell Laboratories in the early 1980s.

The iplots package provides a set of basic plot that are all linked by default.

library(iplots)
with(soil, iplot(northing, resistivity))
with(soil, iplot(easting,northing))

The rggobi package provides another framework for linked plots and brusing within R.

library(rggobi)
ggobi(soil)

Several other tools with linking/brushing support and good integration with R are under development.

• One of the furthest along seems to be rbokeh.
• plotly also provides some support for linked brushing.

## Example: Ethanol Data

The ethanol data frame in the lattice package contains data from an experiment efficiency and emissions in small one-cylinder engines.

The data frame contains 88 observations on three variables:

• NOx: Concentration of nitrogen oxides (NO and NO2) in micrograms.

• C Compression ratio of the engine.

• E Equivalence ratio, a measure of the richness of the air and ethanol fuel mixture.

A scatterplot matrix:

library(lattice)
splom(ethanol)

Conditional plots of NOx against E for each level of C:

xyplot(NOx ~ E | C, data = ethanol)

ggplot(ethanol, aes(x = E, y = NOx)) + geom_point() + facet_wrap(~C)

A single plot with color to separate levels of C:

xyplot(NOx ~ E, group = C, data = ethanol, auto.key = TRUE)

ggplot(ethanol, aes(x = E, y = NOx, color = factor(C))) + geom_point()

Augmented with smooths:

xyplot(NOx ~ E, group = C, type = c("p", "smooth"), data = ethanol,
auto.key = TRUE)

ggplot(ethanol, aes(x = E, y = NOx, color = factor(C))) +
geom_point() + geom_smooth(se = FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Conditional plots of NOx against C given levels of E:

with(ethanol, {
Equivalence.Ratio <- equal.count(E, number = 9, overlap = 0.25)
xyplot(NOx ~ C | Equivalence.Ratio,
panel = function(x, y) {
panel.xyplot(x, y)
panel.loess(x, y, span = 1)
},
aspect = 2.5,
layout = c(5, 2),
xlab = "Compression Ratio",
ylab = "NOx (micrograms/J)")
})

An interactive 3D view:

library(rgl)
clear3d()
with(ethanol,points3d(E, C / 5, NOx / 5,
size = 4, col = rep("black", nrow(ethanol))))

## Example: Earth Quake Locations and Magnitudes

The quakes data frame contains data on locations of seismic events of magnitude 4.0 or larger in a region near Fiji.

The time frame is from 1964 to perhaps 2000. More recent data is available from a number of sources on the web.

Some scatter plot matrices:

library(lattice)
splom(quakes)

splom(quakes, cex = 0.1)

splom(quakes, alpha = 0.1)

splom(quakes, panel = panel.smoothScatter)

Where the quakes are located:

library(maps)
map("world2", c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = "red", cex = 0.1))

map("world2", c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = heat.colors(nrow(quakes))[rank(depth)],
cex = 0.1))

map("world2",c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = heat.colors(nrow(quakes))[rank(depth)]))

Showing depth with grouping:

quakes2 <- mutate(quakes, depth_cut = cut_number(depth, 3))
xyplot(lat ~ long, aspect = "iso", groups = depth_cut,
data = quakes2, auto.key = list(columns = 3, title = "Depth"))

ggplot(quakes2, aes(x = long, y = lat, color = depth_cut) ) +
geom_point() + theme_bw() + coord_map()

ggplot(quakes2, aes(x = long, y = lat, colour = depth_cut, size = mag)) +
geom_point() + theme_bw() + coord_map()

Showing depth with a trellis display:

xyplot(lat ~ long | depth_cut, aspect = "iso", data = quakes2)

xyplot(lat ~ long | depth_cut, aspect = "iso", type=c("p","g"),
data = quakes2, auto.key = TRUE)

ggplot(quakes2, aes(x = long, y = lat)) +
geom_point(color = "blue") + facet_wrap(~depth_cut, nrow = 1) +
theme_bw() + coord_map()

Adding the full data for background context:

xyplot(lat ~ long | cut(depth, 3), aspect = "iso", type=c("p","g"),
panel = function(x, y, ...) {
panel.xyplot(quakes\$long, quakes\$lat, col = "grey")
panel.xyplot(x, y, ...)
},
data = quakes)

## quakes does not contain the depth_cut variable used in the facet
ggplot(quakes2, aes(x = long, y = lat)) +
geom_point(data = quakes, color = "gray") +
geom_point(color = "blue") + facet_wrap(~depth_cut, nrow = 1) +
theme_bw() + coord_map()

3D views of lat, long, and depth:

cloud(-depth ~ long * lat, data = quakes)

library(rgl)
bg3d(col = "black")
clear3d()
with(quakes, points3d(long, lat, -depth / 50, size = 2,
col = rep("white", nrow(quakes))))

clear3d()
with(quakes, points3d(long, lat, -depth / 50, size = 2,
col=heat.colors(nrow(quakes))[rank(mag)]))

## Visualizing 2D and 3D Density Estimates

### 2D Density Estimates

We have already seen some examples of plotting 2D density contours on a scatter plot using geom_density2d.

Here is a version that explicitly creates the density estimate:

data(geyser, package = "MASS")
geyser <- mutate(geyser, duration = lag(duration))
geyser <- filter(geyser, ! is.na(duration))
plot(waiting ~ duration, data = geyser)

sc <- function(x) as.numeric(scale(x))
de <- with(geyser, MASS::kde2d(sc(duration), sc(waiting), n = 50,  h = 0.5))
contour(de)

with(geyser, points(sc(duration), sc(waiting)))

The density estimate can also be viewed with a perspective plot or an interactive 3D plot:

persp(de)

clear3d()
bg3d(col = "white")
surface3d(unique(de\$x), unique(de\$y), de\$z, col = "red")
box3d()
axes3d()

### 3D Density Estimates

A 3D normal mixture:

library(rgl)
library(misc3d)

x <- matrix(rnorm(9000), ncol = 3)
x[1001 : 2000, 2] <- x[1001 : 2000, 2] + 4
x[2001 : 3000, 3] <- x[2001 : 3000, 3] + 4
clear3d()
bg3d(col = "white")
points3d(x = x[,1], y = x[,2], z = x[,3], size = 2, color = "black")