Other 3D Viewing Options

Other options for viewing a 3D scene include:

A stereogram

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)",
            ylab = "Northing (km)", shade=T,
            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)",
            ylab = "Northing (km)", shade=T,
            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)

You must enable Javascript to view this page properly.

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)

You must enable Javascript to view this page properly.

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)

You must enable Javascript to view this page properly.

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)

You must enable Javascript to view this page properly.

clear3d()
spheres3d(dd, radius = 0.05)

You must enable Javascript to view this page properly.

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

You must enable Javascript to view this page properly.

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)))

You must enable Javascript to view this page properly.

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)
## Loading required package: tcltk

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.

A simple example (Rmd).

Linking and Brushing

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:

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))))

You must enable Javascript to view this page properly.

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)
## (loaded the KernSmooth namespace)

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))))

You must enable Javascript to view this page properly.

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

You must enable Javascript to view this page properly.

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()

You must enable Javascript to view this page properly.

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")