plot3D: Plotting multi-dimensional data.

Description Note Note Author(s) References See Also Examples

Description

Functions for visualising 2-D and 3-D data.

Many of the functions are extensions of R's persp or image function.

Other packages that provide visualisation of 3-D data (and which might be better suited) are: rgl, scatterplot3D, misc3D.

Note

This package is dedicated to Carlo.

Note

Some of the functions based on persp will not work properly for all values of phi (which turns the plots upside-down). This is because an assumption is made as to how the perspective plots are viewed.

Author(s)

Karline Soetaert

References

http://www.rforscience.com/rpackages/visualisation/oceanview/

http://www.rforscience.com/rpackages/visualisation/plot3d/

See Also

Functions that are based on the persp function:

Functions defined on the image function:

Other plotting functions:

Colors and colorkey:

Utility functions:

Data sets:

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

Example output

prsp3D> # save plotting parameters
prsp3D>  pm <- par("mfrow")

prsp3D> ## =======================================================================
prsp3D> ## Ribbon, persp, color keys, facets
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D> # simple, no scaling, use breaks to set colors
prsp3D>  persp3D(z = volcano, main = "volcano", clab = c("height", "m"), 
prsp3D+    breaks = seq(80,200, by = 10))

prsp3D> # keep ratios between x- and y (scale = FALSE) 
prsp3D> # change ratio between x- and z (expand)
prsp3D>  persp3D(z = volcano, x = 1: nrow(volcano), y = 1:ncol(volcano),
prsp3D+        expand = 0.3, main = "volcano", facets = FALSE, scale = FALSE,
prsp3D+        clab = "height, m", colkey = list(side = 1, length = 0.5))

prsp3D> # ribbon, in x--direction
prsp3D>  V <- volcano[, seq(1, ncol(volcano), by = 3)]  # lower resolution

prsp3D>  ribbon3D(z = V, colkey = list(width = 0.5, length = 0.5, 
prsp3D+           cex.axis = 0.8, side = 2), clab = "m")

prsp3D> # ribbon, in y-direction
prsp3D>  Vy <- volcano[seq(1, nrow(volcano), by = 3), ]

prsp3D>  ribbon3D(z = Vy, expand = 0.3, space = 0.3, along = "y", 
prsp3D+           colkey = list(width = 0.5, length = 0.5, cex.axis = 0.8))

prsp3D> ## =======================================================================
prsp3D> ## Several ways to visualise 3-D data
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  x <- seq(-pi, pi, by = 0.2)

prsp3D>  y <- seq(-pi, pi, by = 0.3)

prsp3D>  grid <- mesh(x, y)

prsp3D>  z    <- with(grid, cos(x) * sin(y))

prsp3D>  par(mfrow = c(2,2))

prsp3D>  persp3D(z = z, x = x, y = y)

prsp3D>  persp3D(z = z, x = x, y = y, facets = FALSE, curtain = TRUE)

prsp3D> # ribbons in two directions and larger spaces
prsp3D>  ribbon3D(z = z, x = x, y = y, along = "xy", space = 0.3)

prsp3D>  hist3D(z = z, x = x, y = y, border = "black")  

prsp3D> ## =======================================================================
prsp3D> ## Contours and images added
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D>  x <- seq(1, nrow(volcano), by = 3)

prsp3D>  y <- seq(1, ncol(volcano), by = 3) 

prsp3D>  Volcano <- volcano [x, y]

prsp3D>  ribbon3D(z = Volcano, contour = TRUE, zlim= c(-100, 200), 
prsp3D+           image = TRUE)

prsp3D>  persp3D(z = Volcano, contour = TRUE, zlim= c(-200, 200), image = FALSE)

prsp3D>  persp3D(z = Volcano, x = x, y = y, scale = FALSE, 
prsp3D+        contour = list(nlevels = 20, col = "red"), 
prsp3D+        zlim = c(-200, 200), expand = 0.2,
prsp3D+        image = list(col = grey (seq(0, 1, length.out = 100))))

prsp3D>  persp3D(z = Volcano, contour = list(side = c("zmin", "z", "350")), 
prsp3D+        zlim = c(-100, 400), phi = 20, image = list(side = 350))

prsp3D> ## =======================================================================
prsp3D> ## Use of inttype
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D>  persp3D(z = Volcano, shade = 0.5, colkey = FALSE)

prsp3D>  persp3D(z = Volcano, inttype = 2, shade = 0.5, colkey = FALSE)

prsp3D>  x <- y <- seq(0, 2*pi, length.out = 10)

prsp3D>  z <- with (mesh(x, y), cos(x) *sin(y)) + runif(100)

prsp3D>  cv <- matrix(nrow = 10, 0.5*runif(100))

prsp3D>  persp3D(x, y, z, colvar = cv)              # takes averages of z

prsp3D>  persp3D(x, y, z, colvar = cv, inttype = 2) # takes averages of colvar

prsp3D> ## =======================================================================
prsp3D> ## Use of inttype with NAs
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D>  VV <- V2 <- volcano[10:15, 10:15]

prsp3D>  V2[3:4, 3:4] <- NA

prsp3D>  V2[4, 5] <- NA

prsp3D>  image2D(V2, border = "black")  # shows true NA region

prsp3D> # averages of V2, including NAs, NA region larger
prsp3D>  persp3D(z = VV, colvar = V2, inttype = 1, theta = 0, 
prsp3D+        phi = 20, border = "black", main = "inttype = 1")

prsp3D> # extension of VV; NAs unaffected
prsp3D>  persp3D(z = VV, colvar = V2, inttype = 2, theta = 0, 
prsp3D+        phi = 20, border = "black", main = "inttype = 2")

prsp3D> # average of V2, ignoring NA; NA region smaller
prsp3D>  persp3D(z = VV, colvar = V2, inttype = 3, theta = 0, 
prsp3D+        phi = 20, border = "black", main = "inttype = 3")

prsp3D> ## =======================================================================
prsp3D> ## Use of panel.first
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(1, 1))

prsp3D> # A function that is called after the axes were drawn
prsp3D>  panelfirst <- function(trans) {
prsp3D+     zticks <- seq(100, 180, by = 20)
prsp3D+     len <- length(zticks)
prsp3D+     XY0 <- trans3D(x = rep(1, len), y = rep(1, len), z = zticks,
prsp3D+                    pmat = trans)
prsp3D+     XY1 <- trans3D(x = rep(1, len), y = rep(61, len), z = zticks,
prsp3D+                    pmat = trans)
prsp3D+     segments(XY0$x, XY0$y, XY1$x, XY1$y, lty = 2) 
prsp3D+     
prsp3D+     rm <- rowMeans(volcano)             
prsp3D+     XY <- trans3D(x = 1:87, y = rep(ncol(volcano), 87), 
prsp3D+                   z = rm, pmat = trans)
prsp3D+     lines(XY, col = "blue", lwd = 2)
prsp3D+  }

prsp3D>  persp3D(z = volcano, x = 1:87, y = 1: 61, scale = FALSE, theta = 10,
prsp3D+        expand = 0.2, panel.first = panelfirst, colkey = FALSE)

prsp3D> ## =======================================================================
prsp3D> ## with / without colvar / facets 
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D>  persp3D(z = volcano, shade = 0.3, col = gg.col(100))

prsp3D> # shiny colors - set lphi for more brightness
prsp3D>  persp3D(z = volcano, lighting = TRUE, lphi = 90)

prsp3D>  persp3D(z = volcano, col = "lightblue", colvar = NULL, 
prsp3D+    shade = 0.3, bty = "b2")

prsp3D> # this also works:
prsp3D> #  persp3D(z = volcano, col = "grey", shade = 0.3)
prsp3D> 
prsp3D> # tilted x- and y-coordinates of 'volcano'
prsp3D>  volcx <- matrix(nrow = 87, ncol = 61, data = 1:87)

prsp3D>  volcx <- volcx + matrix(nrow = 87, ncol = 61, 
prsp3D+         byrow = TRUE, data = seq(0., 15, length.out = 61))

prsp3D>  volcy <- matrix(ncol = 87, nrow = 61, data = 1:61)

prsp3D>  volcy <- t(volcy + matrix(ncol = 87, nrow = 61, 
prsp3D+         byrow = TRUE, data = seq(0., 15, length.out = 87)))

prsp3D>  persp3D(volcano, x = volcx, y = volcy, phi = 80)

prsp3D> ## =======================================================================
prsp3D> ## Several persps on one plot
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(1, 1))

prsp3D>  clim <- range(volcano)

prsp3D>  persp3D(z = volcano, zlim = c(100, 600), clim = clim, 
prsp3D+    box = FALSE, plot = FALSE)

prsp3D>  persp3D(z = volcano + 200, clim = clim, colvar = volcano, 
prsp3D+        add = TRUE, colkey = FALSE, plot = FALSE)

prsp3D>  persp3D(z = volcano + 400, clim = clim, colvar = volcano, 
prsp3D+        add = TRUE, colkey = FALSE)    # plot = TRUE by default

prsp3D> ## =======================================================================
prsp3D> ## hist3D
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D>  VV   <- volcano[seq(1, 87, 15), seq(1, 61, 15)]

prsp3D>  hist3D(z = VV, scale = FALSE, expand = 0.01, border = "black")

prsp3D> # transparent colors
prsp3D>  hist3D(z = VV, scale = FALSE, expand = 0.01, 
prsp3D+    alpha = 0.5, opaque.top = TRUE, border = "black")

prsp3D>  hist3D(z = VV, scale = FALSE, expand = 0.01, facets = FALSE, lwd = 2)

prsp3D>  hist3D(z = VV, scale = FALSE, expand = 0.01, facets = NA)

prsp3D> ## =======================================================================
prsp3D> ## hist3D and ribbon3D with greyish background, rotated, rescaled,...
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 2))

prsp3D>  hist3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20,
prsp3D+         col = "#0072B2", border = "black", shade = 0.2, ltheta = 90,
prsp3D+         space = 0.3, ticktype = "detailed", d = 2)

prsp3D> # extending the ranges
prsp3D>  plotdev(xlim = c(-0.2, 1.2), ylim = c(-0.2, 1.2), theta = 45)

prsp3D>  ribbon3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20,
prsp3D+         col = "lightblue", border = "black", shade = 0.2, ltheta = 90,
prsp3D+         space = 0.3, ticktype = "detailed", d = 2, curtain = TRUE)

prsp3D>  ribbon3D(z = VV, scale = FALSE, expand = 0.01, bty = "g", phi = 20, zlim = c(95,183),
prsp3D+         col = "lightblue", lighting = TRUE, ltheta = 50, along = "y",
prsp3D+         space = 0.7, ticktype = "detailed", d = 2, curtain = TRUE)

prsp3D> ## =======================================================================
prsp3D> ## hist3D for a 1-D data set
prsp3D> ## =======================================================================
prsp3D> 
prsp3D>  par(mfrow = c(2, 1))

prsp3D>  x <- rchisq(1000, df = 4)

prsp3D>  hs <- hist(x, breaks = 15)

prsp3D>  hist3D(x = hs$mids, y = 1, z = matrix(ncol = 1, data = hs$density), 
prsp3D+   bty = "g", ylim = c(0., 2.0), scale = FALSE, expand = 20, 
prsp3D+   border = "black", col = "white", shade = 0.3, space = 0.1, 
prsp3D+   theta = 20, phi = 20, main = "3-D perspective")

prsp3D> # reset plotting parameters
prsp3D>  par(mfrow = pm)

surf3D> # save plotting parameters
surf3D>  pm   <- par("mfrow")

surf3D>  pmar <- par("mar")

surf3D>  par(mar = c(1, 1, 1, 1))

surf3D> ## =======================================================================
surf3D> ## A three-dimensional shape 
surf3D> ## (ala http://docs.enthought.com/mayavi/mayavi/mlab.html)
surf3D> ## =======================================================================
surf3D> 
surf3D>  par(mfrow = c(2, 2))

surf3D> # create grid matrices
surf3D>  X       <- seq(0, pi, length.out = 50)

surf3D>  Y       <- seq(0, 2*pi, length.out = 50)

surf3D>  M       <- mesh(X, Y)

surf3D>  phi     <- M$x

surf3D>  theta   <- M$y

surf3D> # x, y and z grids
surf3D>  r <- sin(4*phi)^3 + cos(2*phi)^3 + sin(6*theta)^2 + cos(6*theta)^4

surf3D>  x <- r * sin(phi) * cos(theta)

surf3D>  y <- r * cos(phi)

surf3D>  z <- r * sin(phi) * sin(theta)

surf3D> # full colored image
surf3D>  surf3D(x, y, z, colvar = y, colkey = FALSE, shade = 0.5,
surf3D+         box = FALSE, theta = 60)

surf3D> # same, but just facets
surf3D>  surf3D(x, y, z, colvar = y, colkey = FALSE, box = FALSE, 
surf3D+         theta = 60, facets = FALSE)

surf3D> # with colors and border, AND increasing the size
surf3D> # (by reducing the x- y and z- ranges
surf3D>  surf3D(x, y, z, colvar = y, colkey = FALSE, box = FALSE, 
surf3D+         theta = 60, border = "black", xlim = range(x)*0.8, 
surf3D+         ylim = range(y)*0.8, zlim = range(z)*0.8)

surf3D> # Now with one color and shading
surf3D>  surf3D(x, y, z, box = FALSE,
surf3D+         theta = 60, col = "lightblue", shade = 0.9)

surf3D> ## Not run: 
surf3D> ##D  # rotation
surf3D> ##D   for (angle in seq(0, 360, by = 10))
surf3D> ##D     plotdev(theta = angle)
surf3D> ##D 
surf3D> ## End(Not run)
surf3D> 
surf3D> ## =======================================================================
surf3D> ## Several other shapes 
surf3D> ## http://xahlee.info/surface/gallery.html
surf3D> ## =======================================================================
surf3D> 
surf3D>  par(mfrow = c(2, 2)) 

surf3D>  # Shape 1
surf3D>  M  <- mesh(seq(0,  6*pi, length.out = 50), 
surf3D+             seq(pi/3, pi, length.out = 50))

surf3D>  u  <- M$x ; v <- M$y

surf3D>  x <- u/2 * sin(v) * cos(u)

surf3D>  y <- u/2 * sin(v) * sin(u)

surf3D>  z <- u/2 * cos(v)

surf3D>  surf3D(x, y, z, colvar = z, colkey = FALSE, box = FALSE, phi = 50)

surf3D> # Shape 2: add border
surf3D>  M  <- mesh(seq(0, 2*pi, length.out = 50), 
surf3D+             seq(0, 2*pi, length.out = 50))

surf3D>  u  <- M$x ; v  <- M$y

surf3D>  x  <- sin(u)

surf3D>  y  <- sin(v)

surf3D>  z  <- sin(u + v)

surf3D>  surf3D(x, y, z, colvar = z, border = "black", 
surf3D+         colkey = FALSE)

surf3D> # shape 3: uses same mesh, other perspective (d >1)
surf3D>  x <- (3 + cos(v/2)*sin(u) - sin(v/2)*sin(2*u))*cos(v)

surf3D>  y <- (3 + cos(v/2)*sin(u) - sin(v/2)*sin(2*u))*sin(v)

surf3D>  z <- sin(v/2)*sin(u) + cos(v/2)*sin(2*u)

surf3D>  surf3D(x, y, z, colvar = z, colkey = FALSE, d = 2, facets = FALSE)

surf3D> # shape 4: more complex colvar
surf3D>  M  <- mesh(seq(-13.2, 13.2, length.out = 50), 
surf3D+             seq(-37.4, 37.4, length.out = 50))

surf3D>  u  <- M$x   ; v <- M$y

surf3D>  b <- 0.4; r <- 1 - b^2; w <- sqrt(r)

surf3D>  D <- b*((w*cosh(b*u))^2 + (b*sin(w*v))^2)

surf3D>  x <- -u + (2*r*cosh(b*u)*sinh(b*u)) / D

surf3D>  y <- (2*w*cosh(b*u)*(-(w*cos(v)*cos(w*v)) - sin(v)*sin(w*v))) / D

surf3D>  z <- (2*w*cosh(b*u)*(-(w*sin(v)*cos(w*v)) + cos(v)*sin(w*v))) / D

surf3D>  surf3D(x, y, z, colvar = sqrt(x + 8.3), colkey = FALSE, 
surf3D+         theta = 10, border = "black", box = FALSE)

surf3D>  box()

surf3D> ## =======================================================================
surf3D> ## A sphere, with box type with grid lines
surf3D> ## =======================================================================
surf3D> 
surf3D>  par(mar = c(2, 2, 2, 2))

surf3D>  par(mfrow = c(1, 1))

surf3D>  M  <- mesh(seq(0, 2*pi, length.out = 50), 
surf3D+             seq(0,   pi, length.out = 50))

surf3D>  u  <- M$x ; v  <- M$y

surf3D>  x <- cos(u)*sin(v)

surf3D>  y <- sin(u)*sin(v)

surf3D>  z <- cos(v)

surf3D>  colvar <- sin(u*6) * sin(v*6)

surf3D>  surf3D(y, x, z, colvar = colvar, phi = 0, bty = "b2", 
surf3D+         lighting = TRUE, ltheta = 40)

surf3D> ## =======================================================================
surf3D> ## Function spheresurf3D
surf3D> ## =======================================================================
surf3D> 
surf3D>  par(mfrow = c(2, 2))

surf3D>  spheresurf3D()

surf3D> # true ranges are [-1, 1]; set limits to [-0.8, 0.8] to make larger plots
surf3D>  lim <- c(-0.8, 0.8)

surf3D>  spheresurf3D(colkey = FALSE, xlim = lim, ylim = lim, zlim = lim)

surf3D>  spheresurf3D(bty = "b", ticktype = "detailed", phi = 50)

surf3D>  spheresurf3D(colvar = matrix(nrow = 30, data = runif(900)))

surf3D> ## =======================================================================
surf3D> ## Images on a sphere
surf3D> ## =======================================================================
surf3D> 
surf3D>  par(mfrow = c(1, 1), mar = c(1, 1, 1, 3))

surf3D>  AA <- Hypsometry$z; AA[AA<=0] <- NA

surf3D>  lim <- c(-0.8, 0.8)

surf3D> # log transformation of color variable
surf3D>  spheresurf3D(AA, NAcol = "black", theta = 90, phi = 30, box = FALSE,
surf3D+    xlim = lim, ylim = lim, zlim = lim, log = "c")

surf3D> # restore plotting parameters
surf3D>  par(mfrow = pm)

surf3D>  par(mar = pmar)

slic3D> # save plotting parameters
slic3D>  pm <- par("mfrow")

slic3D>  pmar <- par("mar")

slic3D> ## =======================================================================
slic3D> ## Simple slice3D examples
slic3D> ## =======================================================================
slic3D> 
slic3D>  par(mfrow = c(2, 2))

slic3D>  x <- y <- z <- seq(-1, 1, by = 0.1)

slic3D>  grid   <- mesh(x, y, z)

slic3D>  colvar <- with(grid, x*exp(-x^2 - y^2 - z^2))

slic3D> # default is just the panels
slic3D>  slice3D  (x, y, z, colvar = colvar, theta = 60)

slic3D> # contour slices
slic3D>  slicecont3D (x, y, z, ys = seq(-1, 1, by = 0.5), colvar = colvar, 
slic3D+            theta = 60, border = "black")

slic3D>  slice3D  (x, y, z, xs = c(-1, -0.5, 0.5), ys = c(-1, 0, 1), 
slic3D+            zs = c(-1, 0), colvar = colvar, 
slic3D+            theta = 60, phi = 40)

slic3D> ## =======================================================================
slic3D> ## coloring on a surface
slic3D> ## =======================================================================
slic3D> 
slic3D>  XY <- mesh(x, y)

slic3D>  ZZ <- XY$x*XY$y

slic3D>  slice3D  (x, y, z, xs = XY$x, ys = XY$y, zs = ZZ, colvar = colvar, 
slic3D+            lighting =  TRUE, lphi = 90, ltheta = 0)

slic3D> ## =======================================================================
slic3D> ## Specifying transparent colors
slic3D> ## =======================================================================
slic3D> 
slic3D>  par(mfrow = c(1, 1))

slic3D>  x <- y <- z <- seq(-4, 4, by = 0.2)

slic3D>  M <- mesh(x, y, z)

slic3D>  R <- with (M, sqrt(x^2 + y^2 + z^2))

slic3D>  p <- sin(2*R) /(R+1e-3)

slic3D> ## Not run: 
slic3D> ##D # This is very slow - alpha = 0.5 makes it transparent
slic3D> ##D 
slic3D> ##D  slice3D(x, y, z, colvar = p, col = jet.col(alpha = 0.5), 
slic3D> ##D          xs = 0, ys = c(-4, 0, 4), zs = NULL, d = 2) 
slic3D> ## End(Not run)
slic3D> 
slic3D>  slice3D(x, y, z, colvar = p, d = 2, theta = 60, border = "black",
slic3D+          xs = c(-4, 0), ys = c(-4, 0, 4), zs = c(-4, 0))

slic3D> ## =======================================================================
slic3D> ## A section along a transect
slic3D> ## =======================================================================
slic3D> 
slic3D>  data(Oxsat)

slic3D>  Ox <- Oxsat$val[,  Oxsat$lat > - 5 & Oxsat$lat < 5, ]

slic3D>  slice3D(x = Oxsat$lon, z = -Oxsat$depth, y = 1:5, colvar = Ox, 
slic3D+          ys = 1:5, zs = NULL, NAcol = "black", 
slic3D+          expand = 0.4, theta = 45, phi = 45)

slic3D> ## =======================================================================
slic3D> ## isosurf3D example - rather slow
slic3D> ## =======================================================================
slic3D> 
slic3D>  par(mfrow = c(2, 2), mar  = c(2, 2, 2, 2))

slic3D>  x <- y <- z <- seq(-2, 2, length.out = 15)

slic3D>  xyz <- mesh(x, y, z)

slic3D>  F <- with(xyz, log(x^2 + y^2 + z^2 + 
slic3D+                 10*(x^2 + y^2) * (y^2 + z^2) ^2))

slic3D> # use shading for level = 1 - show triangulation with border
slic3D>  isosurf3D(x, y, z, F, level = 1, shade = 0.9, 
slic3D+            col = "yellow", border = "orange")

slic3D> # lighting for level - 2
slic3D>  isosurf3D(x, y, z, F, level = 2, lighting = TRUE,
slic3D+            lphi = 0, ltheta = 0, col = "blue", shade = NA)  

slic3D> # three levels, transparency added
slic3D>  isosurf3D(x, y, z, F, level = seq(0, 4, by = 2), 
slic3D+    col = c("red", "blue", "yellow"), 
slic3D+    clab = "F", alpha = 0.2, theta = 0, lighting = TRUE)  

slic3D> # transparency can also be added afterwards with plotdev()
slic3D> ## Not run: 
slic3D> ##D  isosurf3D(x, y, z, F, level = seq(0, 4, by = 2), 
slic3D> ##D    col = c("red", "blue", "yellow"), 
slic3D> ##D    shade = NA, plot = FALSE, clab = "F")  
slic3D> ##D  plotdev(lighting = TRUE, alpha = 0.2, theta = 0)
slic3D> ## End(Not run)
slic3D> # use of creatisosurf
slic3D>  iso <- createisosurf(x, y, z, F, level = 2)

slic3D>  head(iso)
              x          y         z
[1,]  0.0000000 -0.1179802 -2.000000
[2,] -0.1204879  0.0000000 -2.000000
[3,]  0.0000000 -0.2073773 -1.714286
[4,]  0.0000000 -0.2073773 -1.714286
[5,] -0.1204879  0.0000000 -2.000000
[6,] -0.2138894  0.0000000 -1.714286

slic3D>  triangle3D(iso, col = "green", shade = 0.3)

slic3D> ## Not run: 
slic3D> ##D  # higher resolution
slic3D> ##D   x <- y <- z <- seq(-2, 2, length.out = 50)
slic3D> ##D   xyz <- mesh(x, y, z)
slic3D> ##D   F <- with(xyz, log(x^2 + y^2 + z^2 + 
slic3D> ##D                 10*(x^2 + y^2) * (y^2 + z^2) ^2))
slic3D> ##D 
slic3D> ##D # three levels
slic3D> ##D   isosurf3D(x, y, z, F, level = seq(0, 4, by = 2), 
slic3D> ##D     col = c("red", "blue", "yellow"), 
slic3D> ##D     shade = NA, plot = FALSE, clab = "F")  
slic3D> ##D   plotdev(lighting = TRUE, alpha = 0.2, theta = 0)
slic3D> ## End(Not run)
slic3D> 
slic3D> ## =======================================================================
slic3D> ## voxel3D example
slic3D> ## =======================================================================
slic3D> 
slic3D>  par(mfrow = c(2, 2), mar  = c(2, 2, 2, 2))

slic3D> # fast but needs high resolution grid
slic3D>  x <- y <- z <- seq(-2, 2, length.out = 70)

slic3D>  xyz <- mesh(x, y, z)

slic3D>  F <- with(xyz, log(x^2 + y^2 + z^2 + 
slic3D+                 10*(x^2 + y^2) * (y^2 + z^2) ^2))

slic3D>  voxel3D(x, y, z, F, level = 4, pch = ".", cex = 5)

slic3D> ## =======================================================================
slic3D> ## rotation 
slic3D> ## =======================================================================
slic3D> 
slic3D>  plotdev(theta = 45, phi = 0)

slic3D>  plotdev(theta = 90, phi = 10)

slic3D> # same using createvoxel -  more flexible for coloring
slic3D>  vox <- createvoxel(x, y, z, F, level = 4)

slic3D>  scatter3D(vox$x, vox$y, vox$z, colvar = vox$y, 
slic3D+    bty = "g", colkey = FALSE)

slic3D> ## =======================================================================
slic3D> ## voxel3D to show hypox sites
slic3D> ## =======================================================================
slic3D> 
slic3D>  par(mfrow = c(1, 1), mar = c(2, 2, 2, 2))

slic3D>  Hypox <- createvoxel(Oxsat$lon, Oxsat$lat, Oxsat$depth[1:19], 
slic3D+                       Oxsat$val[,,1:19], level = 40, operator = "<")

slic3D>  panel <- function(pmat) {  # an image at the bottom
slic3D+    Nx <- length(Oxsat$lon)
slic3D+    Ny <- length(Oxsat$lat)
slic3D+    M <- mesh(Oxsat$lon, Oxsat$lat) 
slic3D+    xy <- trans3D(pmat = pmat, x = as.vector(M$x), y = as.vector(M$y), 
slic3D+         z = rep(-1000, length.out = Nx*Ny)) 
slic3D+    x <- matrix(nrow = Nx, ncol = Ny, data = xy$x)
slic3D+    y <- matrix(nrow = Nx, ncol = Ny, data = xy$y)
slic3D+    Bat <- Oxsat$val[,,1]; Bat[!is.na(Bat)] <- 1
slic3D+    image2D(x = x, y = y, z = Bat, NAcol = "black", col = "grey",
slic3D+          add = TRUE, colkey = FALSE)
slic3D+  }

slic3D>  scatter3D(Hypox$x, Hypox$y, -Hypox$z, colvar = Hypox$cv, 
slic3D+            panel.first = panel, pch = ".", bty = "b", 
slic3D+            theta = 30, phi = 20, ticktype = "detailed",
slic3D+            zlim = c(-1000,0), xlim = range(Oxsat$lon), 
slic3D+            ylim = range(Oxsat$lat) )

slic3D> # restore plotting parameters
slic3D>  par(mfrow = pm)

slic3D>  par(mar = pmar)

sctt3D> # save plotting parameters
sctt3D>  pm <- par("mfrow")

sctt3D> ## =======================================================================
sctt3D> ## A sphere 
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 1))

sctt3D>  M  <- mesh(seq(0, 2*pi, length.out = 100), 
sctt3D+             seq(0,   pi, length.out = 100))

sctt3D>  u  <- M$x ; v  <- M$y

sctt3D>  x <- cos(u)*sin(v)

sctt3D>  y <- sin(u)*sin(v)

sctt3D>  z <- cos(v)

sctt3D> # full  panels of box are drawn (bty = "f")
sctt3D>  scatter3D(x, y, z, pch = ".", col = "red", 
sctt3D+            bty = "f", cex = 2, colkey = FALSE)

sctt3D> ## =======================================================================
sctt3D> ## Different types
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par (mfrow = c(2, 2))

sctt3D>  z <- seq(0, 10, 0.2)

sctt3D>  x <- cos(z)

sctt3D>  y <- sin(z)*z

sctt3D> # greyish background for the boxtype (bty = "g") 
sctt3D>  scatter3D(x, y, z, phi = 0, bty = "g",
sctt3D+            pch = 20, cex = 2, ticktype = "detailed")

sctt3D> # add another point
sctt3D>  scatter3D(x = 0, y = 0, z = 0, add = TRUE, colkey = FALSE, 
sctt3D+            pch = 18, cex = 3, col = "black")

sctt3D> # add text
sctt3D>  text3D(x = cos(1:10), y = (sin(1:10)*(1:10) - 1), 
sctt3D+         z = 1:10, colkey = FALSE, add = TRUE, 
sctt3D+         labels = LETTERS[1:10], col = c("black", "red"))

sctt3D> # line plot
sctt3D>  scatter3D(x, y, z, phi = 0, bty = "g", type = "l", 
sctt3D+            ticktype = "detailed", lwd = 4)

sctt3D> # points and lines
sctt3D>  scatter3D(x, y, z, phi = 0, bty = "g", type = "b", 
sctt3D+            ticktype = "detailed", pch = 20, 
sctt3D+            cex = c(0.5, 1, 1.5))

sctt3D> # vertical lines
sctt3D>  scatter3D(x, y, z, phi = 0, bty = "g",  type = "h", 
sctt3D+            ticktype = "detailed")

sctt3D> ## =======================================================================
sctt3D> ## With confidence interval
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  x <- runif(20)

sctt3D>  y <- runif(20)

sctt3D>  z <- runif(20)

sctt3D>  par(mfrow = c(1, 1))

sctt3D>  CI <- list(z = matrix(nrow = length(x),
sctt3D+                        data = rep(0.05, 2*length(x))))

sctt3D> # greyish background for the boxtype (bty = "g")
sctt3D>  scatter3D(x, y, z, phi = 0, bty = "g", CI = CI,
sctt3D+    col = gg.col(100), pch = 18, cex = 2, ticktype = "detailed",
sctt3D+    xlim = c(0, 1), ylim = c(0, 1), zlim = c(0, 1))

sctt3D> # add new set of points
sctt3D>  x <- runif(20)

sctt3D>  y <- runif(20)

sctt3D>  z <- runif(20)

sctt3D>  CI2 <- list(x = matrix(nrow = length(x),
sctt3D+                        data = rep(0.05, 2*length(x))),
sctt3D+              z = matrix(nrow = length(x),
sctt3D+                        data = rep(0.05, 2*length(x))))

sctt3D>  scatter3D(x, y, z, CI = CI2, add = TRUE, col = "red", pch = 16)

sctt3D> ## =======================================================================
sctt3D> ## With a surface
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 1))

sctt3D> # surface = volcano
sctt3D>  M <- mesh(1:nrow(volcano), 1:ncol(volcano))

sctt3D> # 100 points above volcano 
sctt3D>  N  <- 100

sctt3D>  xs <- runif(N) * 87

sctt3D>  ys <- runif(N) * 61

sctt3D>  zs <- runif(N)*50 + 154

sctt3D> # scatter + surface
sctt3D>  scatter3D(xs, ys, zs, ticktype = "detailed", pch = 16, 
sctt3D+    bty = "f", xlim = c(1, 87), ylim = c(1,61), zlim = c(94, 215), 
sctt3D+    surf = list(x = M$x, y = M$y, z = volcano,  
sctt3D+                NAcol = "grey", shade = 0.1))

sctt3D> ## =======================================================================
sctt3D> ## A surface and CI
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 1))

sctt3D>  M <- mesh(seq(0, 2*pi, length = 30), (1:30)/100)

sctt3D>  z <- with (M, sin(x) + y)

sctt3D> # points 'sampled'
sctt3D>  N <- 30

sctt3D>  xs <- runif(N) * 2*pi

sctt3D>  ys <- runif(N) * 0.3

sctt3D>  zs <- sin(xs) + ys + rnorm(N)*0.3

sctt3D>  CI <- list(z = matrix(nrow = length(xs), 
sctt3D+                        data = rep(0.3, 2*length(xs))),
sctt3D+             lwd = 3)

sctt3D> # facets = NA makes a transparent surface; borders are black
sctt3D>  scatter3D(xs, ys, zs, ticktype = "detailed", pch = 16, 
sctt3D+    xlim = c(0, 2*pi), ylim = c(0, 0.3), zlim = c(-1.5, 1.5), 
sctt3D+    CI = CI, theta = 20, phi = 30, cex = 2,
sctt3D+    surf = list(x = M$x, y = M$y, z = z, border = "black", facets = NA)
sctt3D+    )

sctt3D> ## =======================================================================
sctt3D> ## droplines till the fitted surface
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  with (mtcars, {
sctt3D+ 
sctt3D+   # linear regression
sctt3D+    fit <- lm(mpg ~ wt + disp)
sctt3D+ 
sctt3D+   # predict values on regular xy grid
sctt3D+    wt.pred <- seq(1.5, 5.5, length.out = 30)
sctt3D+    disp.pred <- seq(71, 472, length.out = 30)
sctt3D+    xy <- expand.grid(wt = wt.pred, 
sctt3D+                      disp = disp.pred)
sctt3D+ 
sctt3D+    mpg.pred <- matrix (nrow = 30, ncol = 30, 
sctt3D+       data = predict(fit, newdata = data.frame(xy), 
sctt3D+       interval = "prediction"))
sctt3D+ 
sctt3D+ # fitted points for droplines to surface
sctt3D+    fitpoints <- predict(fit) 
sctt3D+ 
sctt3D+    scatter3D(z = mpg, x = wt, y = disp, pch = 18, cex = 2, 
sctt3D+       theta = 20, phi = 20, ticktype = "detailed",
sctt3D+       xlab = "wt", ylab = "disp", zlab = "mpg",  
sctt3D+       surf = list(x = wt.pred, y = disp.pred, z = mpg.pred,  
sctt3D+                   facets = NA, fit = fitpoints),
sctt3D+       main = "mtcars")
sctt3D+   
sctt3D+  })

sctt3D> ## =======================================================================
sctt3D> ## Two ways to make a scatter 3D of quakes data set
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 1)) 

sctt3D> # first way, use vertical spikes (type = "h")
sctt3D>  with(quakes, scatter3D(x = long, y = lat, z = -depth, colvar = mag, 
sctt3D+       pch = 16, cex = 1.5, xlab = "longitude", ylab = "latitude", 
sctt3D+       zlab = "depth, km", clab = c("Richter","Magnitude"),
sctt3D+       main = "Earthquakes off Fiji", ticktype = "detailed", 
sctt3D+       type = "h", theta = 10, d = 2, 
sctt3D+       colkey = list(length = 0.5, width = 0.5, cex.clab = 0.75))
sctt3D+       )

sctt3D> # second way: add dots on bottom and left panel
sctt3D> # before the scatters are drawn, 
sctt3D> # add small dots on basal plane and on the depth plane
sctt3D>  panelfirst <- function(pmat) {
sctt3D+     zmin <- min(-quakes$depth)
sctt3D+     XY <- trans3D(quakes$long, quakes$lat, 
sctt3D+                   z = rep(zmin, nrow(quakes)), pmat = pmat)
sctt3D+     scatter2D(XY$x, XY$y, colvar = quakes$mag, pch = ".", 
sctt3D+             cex = 2, add = TRUE, colkey = FALSE)
sctt3D+  
sctt3D+     xmin <- min(quakes$long)
sctt3D+     XY <- trans3D(x = rep(xmin, nrow(quakes)), y = quakes$lat, 
sctt3D+                   z = -quakes$depth, pmat = pmat)
sctt3D+     scatter2D(XY$x, XY$y, colvar = quakes$mag, pch = ".", 
sctt3D+             cex = 2, add = TRUE, colkey = FALSE)
sctt3D+  }

sctt3D>  with(quakes, scatter3D(x = long, y = lat, z = -depth, colvar = mag, 
sctt3D+       pch = 16, cex = 1.5, xlab = "longitude", ylab = "latitude", 
sctt3D+       zlab = "depth, km", clab = c("Richter","Magnitude"),
sctt3D+       main = "Earthquakes off Fiji", ticktype = "detailed", 
sctt3D+       panel.first = panelfirst, theta = 10, d = 2, 
sctt3D+       colkey = list(length = 0.5, width = 0.5, cex.clab = 0.75))
sctt3D+       )

sctt3D> ## =======================================================================
sctt3D> ## text3D and scatter3D
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  with(USArrests, text3D(Murder, Assault, Rape, 
sctt3D+     colvar = UrbanPop, col = gg.col(100), theta = 60, phi = 20,
sctt3D+     xlab = "Murder", ylab = "Assault", zlab = "Rape", 
sctt3D+     main = "USA arrests", 
sctt3D+     labels = rownames(USArrests), cex = 0.6, 
sctt3D+     bty = "g", ticktype = "detailed", d = 2,
sctt3D+     clab = c("Urban","Pop"), adj = 0.5, font = 2))

sctt3D>  with(USArrests, scatter3D(Murder, Assault, Rape - 1, 
sctt3D+     colvar = UrbanPop, col = gg.col(100), 
sctt3D+     type = "h", pch = ".", add = TRUE))

sctt3D> ## =======================================================================
sctt3D> ## zoom near origin
sctt3D> ## =======================================================================
sctt3D> 
sctt3D> # display axis ranges
sctt3D>  getplist()[c("xlim","ylim","zlim")] 
$xlim
[1]  0.8 17.4

$ylim
[1]  45 337

$zlim
[1]  7.3 46.0


sctt3D> # choose suitable ranges
sctt3D>  plotdev(xlim = c(0, 10), ylim = c(40, 150), 
sctt3D+          zlim = c(7, 25))

sctt3D> ## =======================================================================
sctt3D> ## text3D to label x- and y axis
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 1))

sctt3D>  hist3D (x = 1:5, y = 1:4, z = VADeaths,
sctt3D+         bty = "g", phi = 20,  theta = -60,
sctt3D+         xlab = "", ylab = "", zlab = "", main = "VADeaths",
sctt3D+         col = "#0072B2", border = "black", shade = 0.8,
sctt3D+         ticktype = "detailed", space = 0.15, d = 2, cex.axis = 1e-9)

sctt3D>  text3D(x = 1:5, y = rep(0.5, 5), z = rep(3, 5),
sctt3D+        labels = rownames(VADeaths),
sctt3D+        add = TRUE, adj = 0)

sctt3D>  text3D(x = rep(1, 4),   y = 1:4, z = rep(0, 4),
sctt3D+        labels  = colnames(VADeaths),
sctt3D+        add = TRUE, adj = 1)

sctt3D> ## =======================================================================
sctt3D> ## Scatter2D; bty can also be set = to one of the perspbox alernatives
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(2, 2))

sctt3D>  x <- seq(0, 2*pi, length.out = 30)

sctt3D>  scatter2D(x, sin(x), colvar = cos(x), pch = 16, 
sctt3D+          ylab = "sin", clab = "cos", cex = 1.5)

sctt3D> # other box types:
sctt3D>  scatter2D(x, sin(x), colvar = cos(x), type = "l", lwd = 4, bty = "g")

sctt3D>  scatter2D(x, sin(x), colvar = cos(x), type = "b", lwd = 2, bty = "b2")

sctt3D> # transparent colors and spikes
sctt3D>  scatter2D(x, sin(x), colvar = cos(x), type = "h", lwd = 4, alpha = 0.5)

sctt3D> ## =======================================================================
sctt3D> ## mesh examples and scatter2D
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 2))

sctt3D>  x <- seq(-1, 1, by = 0.1)

sctt3D>  y <- seq(-2, 2, by = 0.2)

sctt3D>  grid <- mesh(x, y)

sctt3D>  z    <- with(grid, cos(x) * sin(y))

sctt3D>  image2D(z, x = x, y = y)

sctt3D>  points(grid)

sctt3D>  scatter2D(grid$x, grid$y, colvar = z, pch = 20, cex = 2)

sctt3D> ## =======================================================================
sctt3D> ## scatter plot with confidence intervals
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(2, 2))

sctt3D>  x  <- sort(rnorm(10)) 

sctt3D>  y  <- runif(10)

sctt3D>  cv <- sqrt(x^2 + y^2)

sctt3D>  CI <- list(lwd = 2)

sctt3D>  CI$x <- matrix (nrow = length(x), data = c(rep(0.25, 2*length(x))))

sctt3D>  scatter2D(x, y, colvar = cv, pch = 16, cex = 2, CI = CI)

sctt3D>  scatter2D(x, y, colvar = cv, pch = 16, cex = 2, CI = CI, type = "b")

sctt3D>  CI$y <- matrix (nrow = length(x), data = c(rep(0.05, 2*length(x))))

sctt3D>  CI$col <- "black"

sctt3D>  scatter2D(x, y, colvar = cv, pch = 16, cex = 2, CI = CI)

sctt3D>  CI$y[c(2,4,8,10), ] <- NA  # Some points have no CI

sctt3D>  CI$x[c(2,4,8,10), ] <- NA  # Some points have no CI

sctt3D>  CI$alen <- 0.02            # increase arrow head

sctt3D>  scatter2D(x, y, colvar = cv, pch = 16, cex = 2, CI = CI)

sctt3D> ## =======================================================================
sctt3D> ## Scatter on an image
sctt3D> ## =======================================================================
sctt3D>  
sctt3D>  par(mfrow = c(1, 1))

sctt3D> # image of oxygen saturation
sctt3D>  oxlim <- range(Oxsat$val[,,1], na.rm  = TRUE) 

sctt3D>  image2D(z = Oxsat$val[,,1], x = Oxsat$lon, y = Oxsat$lat,
sctt3D+        contour = TRUE, 
sctt3D+        xlab = "longitude", ylab = "latitude", 
sctt3D+        main = "Oxygen saturation", clim = oxlim, clab = "%")

sctt3D> # (imaginary) measurements at 5 sites
sctt3D>  lon   <- c( 11.2,   6.0, 0.9,  -4, -8.8)

sctt3D>  lat   <- c(-19.7,-14.45,-9.1,-3.8, -1.5)

sctt3D>  O2sat <- c(   90,    95,  92,  85,  100)

sctt3D> # add to image; use same zrange; avoid adding  a color key
sctt3D>  scatter2D(colvar = O2sat, x = lon, y = lat, clim = oxlim, pch = 16,
sctt3D+          add = TRUE, cex = 2, colkey = FALSE)

sctt3D> ## =======================================================================
sctt3D> ## Scatter on a contourplot
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  par(mfrow = c(1, 1))

sctt3D> # room for colorkey by setting colkey = list(plot = FALSE)
sctt3D> 
sctt3D> # contour plot of the ocean's bathymetry
sctt3D>  Depth <- Hypsometry$z

sctt3D>  Depth[Depth > 0] <- NA

sctt3D>  contour2D(z = Depth, x = Hypsometry$x, y = Hypsometry$y, 
sctt3D+        xlab = "longitude", ylab = "latitude", 
sctt3D+        col = "black", NAcol = "grey", levels = seq(-6000, 0, by = 2000),
sctt3D+        main = "Oxygen saturation along ship track", 
sctt3D+        colkey = list(plot = FALSE))

sctt3D> # add data to image; with  a color key
sctt3D>  scatter2D(colvar = O2sat, x = lon, y = lat, pch = 16,
sctt3D+          add = TRUE, cex = 2, clab = "%")

sctt3D> ## =======================================================================
sctt3D> ## scatter2D for time-series plots
sctt3D> ## =======================================================================
sctt3D> 
sctt3D> # Plotting sunspot 'anomalies'
sctt3D> sunspot <- data.frame(year = time(sunspot.month), 
sctt3D+   anom = sunspot.month - mean(sunspot.month))

sctt3D> # long-term moving average of anomaly
sctt3D> ff <- 100

sctt3D> sunspot$ma <- filter(sunspot$anom, rep(1/ff, ff), sides = 2)

sctt3D> with (sunspot, lines2D(year, anom, 
sctt3D+   colvar = anom > 0, 
sctt3D+   col = c("pink", "lightblue"),
sctt3D+   main = "sunspot anomaly", type = "h", 
sctt3D+   colkey = FALSE, las = 1, xlab = "year", ylab = ""))

sctt3D> lines2D(sunspot$year, sunspot$ma, add = TRUE)  

sctt3D> # The same
sctt3D> #with (sunspot, plot(year, anom, 
sctt3D> #  col = c("pink", "lightblue")[(anom > 0) + 1],
sctt3D> #  main = "sunspot", type = "h", las = 1))
sctt3D> 
sctt3D> # but this does not work due to NAs...  
sctt3D> # lines(sunspot$year, sunspot$ma)  
sctt3D> 
sctt3D> ## =======================================================================
sctt3D> ## text2D
sctt3D> ## =======================================================================
sctt3D> 
sctt3D>  with(USArrests, text2D(x = Murder, y = Assault + 5, colvar = Rape, 
sctt3D+      xlab = "Murder", ylab = "Assault", clab = "Rape", 
sctt3D+      main = "USA arrests", labels = rownames(USArrests), cex = 0.6, 
sctt3D+      adj = 0.5, font = 2))

sctt3D>  with(USArrests, scatter2D(x = Murder, y = Assault, colvar = Rape, 
sctt3D+      pch = 16, add = TRUE, colkey = FALSE))

sctt3D> # reset plotting parameters
sctt3D>  par(mfrow = pm)
Warning message:
In cbind(Col[-1], Col[-len]) :
  number of rows of result is not a multiple of vector length (arg 1)

sgmn3D> # save plotting parameters
sgmn3D>   pm <- par("mfrow")

sgmn3D> ## ========================================================================
sgmn3D> ## arrows, points, segments, box
sgmn3D> ## ========================================================================
sgmn3D> 
sgmn3D> # Create a grid of x, y, and z values
sgmn3D>  xx <- yy <- seq(-0.8, 0.8, by = 0.2)

sgmn3D>  zz <- seq(-0.8, 0.8, by = 0.8)

sgmn3D>  M <- mesh(xx, yy, zz)

sgmn3D>  x0 <- M$x; y0 <- M$y; z0 <- M$z

sgmn3D>  x1 <- x0 + 0.1

sgmn3D>  Col <- c("red", "blue", "green") 

sgmn3D>  arrows3D(x0, y0, z0, x1 = x1, colvar = z0, lwd = 2, 
sgmn3D+           d = 2, clab = "z-value", col = Col, length = 0.1,
sgmn3D+           xlim = c(-0.8, 0.8), ylim = c(-0.8, 0.8),
sgmn3D+           main = "arrows3D, points3D, segments3D, border3D")

sgmn3D> # add starting point of arrows
sgmn3D>  points3D(x0, y0, z0, add = TRUE, colvar = z0, 
sgmn3D+            colkey = FALSE, pch = ".", cex = 3, col = Col)

sgmn3D> # use segments to add section
sgmn3D>  x0 <- c(-0.8, 0.8,  0.8, -0.8)

sgmn3D>  x1 <- c( 0.8, 0.8, -0.8, -0.8)

sgmn3D>  y0 <- c(-0.8, -0.8, 0.8, -0.8)

sgmn3D>  y1 <- c(-0.8,  0.8, 0.8, 0.8)

sgmn3D>  z0 <- c(0., 0., 0., 0.)

sgmn3D>  segments3D(x0, y0, z0, x1, y1, z1 = z0,
sgmn3D+      add = TRUE, col = "black", lwd = 2)

sgmn3D> # add a box 
sgmn3D>  border3D(-0.8, -0.8, -0.8, 0.8, 0.8, 0.8,
sgmn3D+        col = "orange", add = TRUE, lwd = 3)

sgmn3D> ## ========================================================================
sgmn3D> ## boxes, cubes
sgmn3D> ## ========================================================================
sgmn3D> 
sgmn3D> # borders are boxes without facets  
sgmn3D>  border3D(x0 = seq(-0.8, -0.1, by = 0.1), 
sgmn3D+        y0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+        z0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+        x1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+        y1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+        z1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+        col = gg.col(8), lty = 2, 
sgmn3D+        lwd = c(1, 4), phi = 20, main = "border3D")

sgmn3D>  box3D(x0 = -0.8, y0 = -0.8, z0 = -0.8, 
sgmn3D+        x1 = 0.8, y1 = 0.8, z1 = 0.8, 
sgmn3D+        border = "black", lwd = 2, 
sgmn3D+        col = gg.col(1, alpha = 0.8), 
sgmn3D+        main = "box3D")

sgmn3D>  box3D(x0 = seq(-0.8, -0.1, by = 0.1), 
sgmn3D+        y0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+        z0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+        x1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+        y1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+        z1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+        col = rainbow(n = 8, alpha = 0.1), 
sgmn3D+        border = "black", lwd = 2, phi = 20)

sgmn3D> # here the perspective does not always work 
sgmn3D> # use alpha.col to set the transparency of a vector of colors
sgmn3D>  box3D(x0 = runif(3), y0 = runif(3), z0 = runif(3),
sgmn3D+        x1 = runif(3), y1 = runif(3), z1 = runif(3),
sgmn3D+        col = c("red", "lightblue", "orange"), alpha = 0.5,
sgmn3D+        border = "black", lwd = 2)

sgmn3D> ## ========================================================================
sgmn3D> ## rectangles
sgmn3D> ## ========================================================================
sgmn3D> # at constant 'z'
sgmn3D>  rect3D(x0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+         y0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+         z0 = seq(-0.8, -0.1, by = 0.1),
sgmn3D+         x1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+         y1 = seq(0.8, 0.1, by = -0.1),
sgmn3D+         col = gg.col(8), border = "black",
sgmn3D+         bty = "g", lwd = 2, phi = 20, main = "rect3D")

sgmn3D> # constant y and with transparent facets
sgmn3D>  rect3D(x0 = 0, y0 = 0, z0 = 0, x1 = 1, z1 = 5,
sgmn3D+         ylim = c(0, 1), facets = NA, border = "red",
sgmn3D+         bty = "g", lwd = 2, phi = 20)

sgmn3D> # add rect at constant z, with colored facet
sgmn3D>  rect3D(x0 = 0, y0 = 0, z0 = 0, x1 = 1, y1 = 1,
sgmn3D+         border = "red", add = TRUE)

sgmn3D> ## ========================================================================
sgmn3D> ## arrows added to a persp plot 
sgmn3D> ## ========================================================================
sgmn3D> 
sgmn3D>  x <- y <- seq(-10, 10, length = 30)

sgmn3D>  z <- outer(x, y, FUN = function(x,y) x^2 + y^2)

sgmn3D>  persp3D(x, y, z, theta = 30, phi = 30, 
sgmn3D+          col = "lightblue", ltheta = 120, shade = 0.75, 
sgmn3D+          ticktype = "detailed", xlab = "X", 
sgmn3D+          ylab = "Y", zlab = "x^2+y^2" )  

sgmn3D> # Points where to put the arrows
sgmn3D>  x <- y <- seq(-10, 10, len = 6)

sgmn3D>  X0 <- outer(x, y, FUN = function (x,y) x)

sgmn3D>  Y0 <- outer(x, y, FUN = function (x,y) y)

sgmn3D>  Z0 <- outer(x, y, FUN = function (x,y) x^2 + y^2)

sgmn3D>  X1 <- X0 + 1

sgmn3D>  Y1 <- Y0 + 1

sgmn3D>  Z1 <- Z0 + 10

sgmn3D>  arrows3D(X0, Y0, Z0, X1, Y1, Z1, lwd = 2, 
sgmn3D+          add = TRUE, type = "curved", col = "red")

sgmn3D>  segments3D(X0, Y0, Z0, X0, Y0, rep(0, length(X0)), lwd = 2, 
sgmn3D+          add = TRUE, col = "green")

sgmn3D> ## ========================================================================
sgmn3D> ## polygon3D 
sgmn3D> ## ========================================================================
sgmn3D> 
sgmn3D>  x <- runif(10)

sgmn3D>  y <- runif(10)

sgmn3D>  z <- runif(10)

sgmn3D>  polygon3D(x, y, z)

sgmn3D> # several polygons, separated by NAs
sgmn3D>  x <- runif(39) 

sgmn3D>  y <- runif(39)

sgmn3D>  z <- runif(39)

sgmn3D>  ii <- seq(4, 36, by  = 4)

sgmn3D>  x[ii] <- y[ii] <- z[ii] <- NA 

sgmn3D> # transparent colors (alpha)
sgmn3D>  polygon3D(x, y, z, border = "black", lwd = 3,
sgmn3D+    col = gg.col(length(ii) + 1, alpha = 0.8), 
sgmn3D+    main = "polygon3D")

sgmn3D> ## ========================================================================
sgmn3D> ## 2D examples, with color key
sgmn3D> ## ========================================================================
sgmn3D> 
sgmn3D> arrows2D(x0 = runif(10), y0 = runif(10),
sgmn3D+          x1 = runif(10), y1 = runif(10), colvar = 1:10, 
sgmn3D+          code = 3, main = "arrows2D, segments2D")

sgmn3D> segments2D(x0 = runif(10), y0 = runif(10),
sgmn3D+          x1 = runif(10), y1 = runif(10), colvar = 1:10, 
sgmn3D+          lwd = 2, add = TRUE, colkey = FALSE) 

sgmn3D> # transparency
sgmn3D> rect2D(x0 = runif(10), y0 = runif(10),
sgmn3D+        x1 = runif(10), y1 = runif(10), colvar = 1:10, 
sgmn3D+        alpha = 0.4, lwd = 2, main = "rect2D") 

sgmn3D> ## ========================================================================
sgmn3D> ## polygon2D 
sgmn3D> ## ========================================================================
sgmn3D> 
sgmn3D>  x <- runif(10)

sgmn3D>  y <- runif(10)

sgmn3D>  polygon2D(x, y)    # same as polygon

sgmn3D> # several polygons, separated by NAs
sgmn3D>  x <- runif(59) 

sgmn3D>  y <- runif(59)

sgmn3D>  ii <- seq(5, 55, by  = 5)

sgmn3D>  x[ii] <- y[ii] <- NA 

sgmn3D> # transparent colors (alpha)
sgmn3D>  polygon2D(x, y, border = "black", lwd = 3,
sgmn3D+    colvar = 1:(length(ii) + 1), 
sgmn3D+    col = gg.col(), alpha = 0.2,
sgmn3D+    main = "polygon2D")

sgmn3D> # restore plotting parameters
sgmn3D>  par(mfrow = pm)

imag2D> # save plotting parameters
imag2D>  pm <- par("mfrow")

imag2D> ## =======================================================================
imag2D> ## Difference between x or y a vector/matrix and rasterImage
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(2, 2))

imag2D>  x <- y <- 1:3

imag2D>  z <- matrix (nrow = 3, ncol = 3, data = 1:9)

imag2D>  image2D(z, x, y, border = "black")

imag2D>  image2D(z, x, y, rasterImage = TRUE, border = "black")

imag2D>  image2D(z, x = matrix(nrow = 3, ncol = 3, data = x), 
imag2D+        y, border = "black")

imag2D>  image2D(z, x, y, border = "black", theta = 45)

imag2D> ## =======================================================================
imag2D> ## shading, light, adding contours, points and lines
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(2, 2))

imag2D>  nr <- nrow(volcano)

imag2D>  nc <- ncol(volcano)

imag2D>  image2D(volcano, x = 1:nr, y = 1:nc, lighting = TRUE, 
imag2D+        main = "volcano", clab = "height, m")

imag2D>  abline(v = seq(10, 80, by = 10))

imag2D>  abline(h = seq(10, 60, by = 10))

imag2D>  points(50, 30, pch = 3, cex = 5, lwd = 3, col = "white")

imag2D>  image2D(z = volcano, x = 1:nr, y = 1:nc, lwd = 2, shade = 0.2,
imag2D+        main = "volcano", clab = "height, m")

imag2D>  image2D(volcano, x = 1:nr, y = 1:nc, contour = TRUE, shade = 0.5, lphi = 0,
imag2D+        col = "lightblue", main = "volcano")

imag2D>  breaks <- seq(90, 200, by = 10)

imag2D>  image2D(volcano, x = 1:nr, y = 1:nc, col = jet.col(length(breaks)-1),
imag2D+        main = "volcano", clab = "height, m", breaks = breaks)

imag2D> ## =======================================================================
imag2D> ## Contour plots
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(2, 2))

imag2D>  V <- volcano - 150

imag2D> # default, no color key
imag2D>  contour2D(z = V, colkey = FALSE, lwd = 2)

imag2D> # imposed levels
imag2D>  contour2D(z = V, lwd = 2, levels = seq(-40, 40, by = 20))

imag2D> # negative levels dashed
imag2D>  contour2D(z = V, col = "black", lwd = 2, 
imag2D+            levels = seq(0, 40, by = 20))

imag2D>  contour2D(z = V, col = "black", lwd = 2, lty = 2, 
imag2D+            levels = seq(-40, -20, by = 20), add = TRUE)

imag2D> # no labels, imposed number of levels, colorkey
imag2D>  contour2D(z = V, lwd = 2, nlevels = 20, drawlabels = FALSE, 
imag2D+            colkey = list(at = seq(-40, 40, by = 20)))

imag2D> ## =======================================================================
imag2D> ## A large data set, input is an array
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(1, 1))

imag2D>  image2D(z = Oxsat$val[, , 1], x = Oxsat$lon, y = Oxsat$lat,
imag2D+        main = "surface oxygen saturation data 2005", NAcol = "black", 
imag2D+        clab = c("","","%"))

imag2D> # images at first 9 depths - use subset to select them
imag2D>  image2D(z = Oxsat$val, subset = 1:9, 
imag2D+        x = Oxsat$lon, y = Oxsat$lat,
imag2D+        margin = c(1, 2), NAcol = "black", 
imag2D+        xlab = "longitude", ylab = "latitude", 
imag2D+        zlim = c(0, 115),
imag2D+        main = paste("depth ", Oxsat$depth[1:9], " m"),
imag2D+        mfrow = c(3, 3))

imag2D> # images at latitude - depth section - increase resolution
imag2D>  z <- Oxsat$val[,  Oxsat$lat > - 5 & Oxsat$lat < 5, ]

imag2D>  image2D(z = z, x = Oxsat$lon, y = Oxsat$depth,
imag2D+        margin = c(1, 3), NAcol = "black", 
imag2D+        resfac = 3, ylim = c(5000, 0))

imag2D> # show position of transects 
imag2D>  image2D(z = Oxsat$val[ , ,1], 
imag2D+        x = Oxsat$lon, y = Oxsat$lat,
imag2D+        NAcol = "black")

imag2D>  abline(h = Oxsat$lat[Oxsat$lat > - 5 & Oxsat$lat < 5])      

imag2D> ## =======================================================================
imag2D> ## Image of a list of matrices
imag2D> ## =======================================================================
imag2D> 
imag2D>  listvolcano <- list(volcano = volcano, logvolcano = log(volcano)) 

imag2D>  image2D(listvolcano, x = 1:nr, y = 1:nc, contour = TRUE,
imag2D+        main = c("volcano", "log(volcano)"), 
imag2D+        clab = list("height, m", "log(m)"),
imag2D+        zlim = list(c(80, 200), c(4.4, 5.5)))

imag2D> ## =======================================================================
imag2D> ## Image of a list of arrays
imag2D> ## =======================================================================
imag2D> 
imag2D> ## Not run: 
imag2D> ##D # crude conversion from oxsat to oxygen 
imag2D> ##D  listoxygen <- list(Oxsat$val, Oxsat$val/100 * 360)
imag2D> ##D   
imag2D> ##D  image2D(z = listoxygen, 
imag2D> ##D        x = Oxsat$lon, y = Oxsat$lat,
imag2D> ##D        margin = c(1, 2), NAcol = "black", 
imag2D> ##D        main = c("Oxygen saturation ", " Oxygen concentration"),
imag2D> ##D        mtext = paste("depth ", Oxsat$depth, " m")
imag2D> ##D        )
imag2D> ## End(Not run)
imag2D> 
imag2D> ## =======================================================================
imag2D> ## 'x', 'y' and 'z' are matrices
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(2, 1))

imag2D> # tilted x- and y-coordinates of 'volcano'
imag2D>  volcx <- matrix(nrow = 87, ncol = 61, data = 1:87)

imag2D>  volcx <- volcx + matrix(nrow = 87, ncol = 61, 
imag2D+         byrow = TRUE, data = seq(0., 15, length.out = 61))

imag2D>  volcy <- matrix(ncol = 87, nrow = 61, data = 1:61)

imag2D>  volcy <- t(volcy + matrix(ncol = 87, nrow = 61, 
imag2D+         byrow = TRUE, data = seq(0., 25, length.out = 87)))

imag2D>  image2D(volcano, x = volcx, y = volcy)

imag2D> # x and y can also be of dimension dim(z)+1:
imag2D> ## Not run: 
imag2D> ##D # tilted x- and y-coordinates of 'volcano'
imag2D> ##D  volcx <- matrix(nrow = 88, ncol = 62, data = 1:88)
imag2D> ##D  volcx <- volcx + matrix(nrow = 88, ncol = 62,
imag2D> ##D         byrow = TRUE, data = seq(0., 15, length.out = 62))
imag2D> ##D 
imag2D> ##D  volcy <- matrix(ncol = 88, nrow = 62, data = 1:62)
imag2D> ##D  volcy <- t(volcy + matrix(ncol = 88, nrow = 62,
imag2D> ##D         byrow = TRUE, data = seq(0., 25, length.out = 88)))
imag2D> ##D 
imag2D> ##D  image2D(volcano, x = volcx, y = volcy)
imag2D> ## End(Not run)
imag2D> 
imag2D> # use of panel function
imag2D>  image2D(volcano, x = volcx, y = volcy, NAcol = "black", 
imag2D+        panel.first = substitute(box(col = "lightgrey", lwd = 30)))

imag2D> ## =======================================================================
imag2D> ## Image with NAs and logs
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(2, 2))

imag2D> # normal volcano
imag2D>  image2D(volcano, clab = c("height", "m"))

imag2D> # logarithmic z-axis
imag2D>  image2D(volcano, log = "z", clab = c("height", "m"),
imag2D+      main = "log='z'")

imag2D> # Including NAs
imag2D>  VOLC <- volcano - 110

imag2D>  VOLC [VOLC <= 0] <- NA

imag2D>  image2D(VOLC, main = "including NAs and rescaled")

imag2D> # both
imag2D>  image2D(VOLC, NAcol = "black", log = "z", zlim = c(1, 100),
imag2D+      main = "NAs and log = 'z'")

imag2D> ## =======================================================================
imag2D> ## Image with contour specification (alpha sets the transparency)
imag2D> ## =======================================================================
imag2D> 
imag2D>  par(mfrow = c(1, 1))

imag2D>  image2D(volcano, shade = 0.2, rasterImage = TRUE,
imag2D+    contour = list(col = "white", labcex = 0.8, lwd = 3, alpha = 0.5))

imag2D> # same:
imag2D> ## Not run: 
imag2D> ##D  image2D(z = volcano, shade = 0.2, rasterImage = TRUE)
imag2D> ##D  contour2D(z = volcano, col = "white", labcex = 0.8, 
imag2D> ##D    lwd = 3, alpha = 0.5, add = TRUE)
imag2D> ## End(Not run)
imag2D> # reset plotting parameters
imag2D>  par(mfrow = pm)

imag3D> # save plotting parameters
imag3D>  pm <- par("mfrow")

imag3D> ## =======================================================================
imag3D> ## images in x, y, z plane
imag3D> ## =======================================================================
imag3D> 
imag3D>  par(mfrow = c(2, 2))

imag3D> # images in x, y, z plane
imag3D> # We use colkey = list(plot = FALSE) to create room for a color key
imag3D>  image3D(y = seq(0, 1, 0.1), z = seq(0, 1, 0.1), x = 0.5, 
imag3D+    col = "blue", xlim = c(0,1), colkey = list(plot = FALSE))

imag3D>  image3D(x = seq(0, 1, 0.1), z = seq(0, 1, 0.1), y = 0.5, 
imag3D+    add = TRUE, col = "red", alpha = 0.2)   # alpha makes it transparent

imag3D>  image3D(x = seq(0, 1, 0.1), y = seq(0, 1, 0.1), z = 0.5, 
imag3D+    add = TRUE, col = "green")

imag3D>  colkey(col = c("green", "red", "blue"), clim = c(0.5, 3.5), 
imag3D+    at = 1:3, labels = c("z", "y", "x"), add = TRUE)

imag3D> #
imag3D>  image3D(z = 100, colvar = volcano, zlim = c(0, 150),
imag3D+    clab = c("height", "m"))

imag3D> #
imag3D>  image3D( x = 0.5, colvar = volcano, xlim = c(0, 1), 
imag3D+    ylim = c(0, 1), zlim = c(0, 1))

imag3D>  image3D( y = 0.5, colvar = volcano, add = TRUE)

imag3D> #
imag3D>  image3D( z = 1, colvar = volcano, 
imag3D+    x = seq(0, 1, length.out = nrow(volcano)),
imag3D+    y = seq(0, 1, length.out = ncol(volcano)), 
imag3D+    xlim = c(0, 2), ylim = c(0, 2), zlim = c(0, 2))

imag3D>  image3D(y = 2, colvar = volcano, add = TRUE, 
imag3D+     shade = 0.2,
imag3D+     x = seq(0, 1, length.out = nrow(volcano)),
imag3D+     z = seq(1, 2, length.out = ncol(volcano)))

imag3D>  image3D(x = 2, colvar = NULL, col = "orange", add = TRUE, 
imag3D+     y = seq(0, 1, length.out = nrow(volcano)),
imag3D+     z = seq(1, 2, length.out = ncol(volcano)))

imag3D> # reset plotting parameters
imag3D>  par(mfrow = pm)

cntr3D> # save plotting parameters
cntr3D>  pm <- par("mfrow")

cntr3D> ## =======================================================================
cntr3D> ## Contours
cntr3D> ## =======================================================================
cntr3D>  par (mfrow = c(2, 2))

cntr3D>  r <- 1:nrow(volcano)

cntr3D>  c <- 1:ncol(volcano)

cntr3D>  contour3D(x = r, y = c, z = 100, colvar = volcano, zlim = c(0, 150),
cntr3D+    clab = c("height", "m"))

cntr3D>  contour3D(x = 100, y = r, z = c, colvar = volcano, clab = c("height", "m"))

cntr3D>  contour3D(z = volcano, colvar = volcano, lwd = 2, 
cntr3D+    nlevels = 20, clab = c("height", "m"), colkey = FALSE)

cntr3D>  contour3D(y = volcano, colvar = volcano, lwd = 2, 
cntr3D+    nlevels = 10, clab = c("height", "m"))

cntr3D> ## =======================================================================
cntr3D> ## Composite images and contours in 3D
cntr3D> ## =======================================================================
cntr3D>  persp3D(z = volcano, zlim = c(90, 300), col = "white", 
cntr3D+          shade = 0.1, d = 2, plot = FALSE)

cntr3D>  contour3D(z = volcano, colvar = volcano, lwd = 2, add = TRUE, 
cntr3D+          nlevels = 20, clab = c("height", "m"), plot = FALSE,
cntr3D+          colkey = list(at = seq(90, 190, length.out = 5)))

cntr3D>  contour3D(z = 300, colvar = volcano, lwd = 2, col = "grey",
cntr3D+          add = TRUE, nlevels = 5)

cntr3D> ## =======================================================================
cntr3D> ## the viewing depth of contours (dDepth)
cntr3D> ## =======================================================================
cntr3D> 
cntr3D> # too low
cntr3D>  persp3D(z = volcano, col = "white", shade = 0.1, plot = FALSE)

cntr3D>  contour3D(z = volcano, colvar = volcano, lwd = 2, 
cntr3D+          add = TRUE, dDepth = 0, col = "black")

cntr3D> # default
cntr3D>  persp3D(z = volcano, col = "white", shade = 0.1, plot = FALSE)

cntr3D>  contour3D(z = volcano, colvar = volcano, lwd = 2, 
cntr3D+          add = TRUE, dDepth = 0.1, col = "black")

cntr3D> # too high
cntr3D>  persp3D(z = volcano, col = "white", shade = 0.1, plot = FALSE)

cntr3D>  contour3D(z = volcano, colvar = volcano, lwd = 1, 
cntr3D+          add = TRUE, dDepth = 0.5, col = "black")

cntr3D> # reset plotting parameters
cntr3D>  par(mfrow = pm)

colkey> # save plotting parameters
colkey>  pm <- par(mfrow = c(2, 2))

colkey>  pmar <- par(mar = c(5.1, 4.1, 4.1, 2.1))

colkey> ## =======================================================================
colkey> ##  colorkey as argument of a plot3D function
colkey> ## =======================================================================
colkey> # default, colkey = NULL: adds colkey because multiple colors
colkey>  image2D(z = volcano)  

colkey> # default, colkey = NULL: no colkey because only one color
colkey>  image2D(z = volcano, col = "grey", shade = 0.2, contour = TRUE)  

colkey> # colkey = FALSE: no color key, no extra space foreseen
colkey>  image2D(z = volcano, colkey = FALSE)

colkey> # colkey = list(plot = FALSE): no color key, extra space foreseen
colkey>  image2D(z = volcano, colkey = list(plot = FALSE, side = 3))

colkey>  colkey (side = 3, add = TRUE, clim = range(volcano))

colkey> ## =======================================================================
colkey> ##  colorkey in new plot
colkey> ## =======================================================================
colkey>  
colkey>  colkey(side = 1, clim = c(0, 1), add = FALSE, clab = "z", 
colkey+    col.clab = "red", adj.clab = 0)

colkey>  colkey(side = 2, clim = c(0, 1), clab = "z", length = 0.5, width = 0.5)

colkey>  colkey(side = 3, clim = c(0, 1), lwd = 3, clab = c("a","b","c","d"), 
colkey+    line.clab = 5)

colkey>  colkey(side = 4, clim = c(1e-6, 1), clog = TRUE, 
colkey+    clab = "a very long title in bold and close to the key", 
colkey+    line.clab = 1, side.clab = 2, font.clab = 2)

colkey> ## =======================================================================
colkey> ##  colorkey added to existing plot
colkey> ## =======================================================================
colkey> 
colkey>  par(mfrow = c(1, 1))

colkey>  image2D(volcano, xlab = "", clab = "m", 
colkey+        colkey = list(side = 1, length = 0.5, width = 0.5, 
colkey+          line.clab = 1))

colkey>  colkey(side = 3, clim = range(volcano), add = TRUE)

colkey> # 'dist' to put colkey within the image 
colkey> # 'shift' to position colkey to the right or upward
colkey>  par(mfrow = c(1, 1))

colkey>  image2D(volcano, colkey = FALSE)

colkey>  colkey(clim = range(volcano), dist = -0.15, shift = 0.2,
colkey+         side = 3, add = TRUE, clab = "key 1", col.clab = "white",
colkey+         length = 0.5, width = 0.5, col.axis = "white", 
colkey+         col.ticks = "white", cex.axis = 0.8)

colkey>  colkey(clim = range(volcano), dist = -0.1, shift = -0.2,
colkey+         side = 4, add = TRUE, clab = "key 2", col.clab = "white",
colkey+         length = 0.3, width = 0.5, col.axis = "white", 
colkey+         col.ticks = "white", col.box = "red", cex.axis = 0.8)

colkey>  colkey(clim = range(volcano), dist = -0.3, 
colkey+         side = 1, add = TRUE, clab = "key 3", col.clab = "white",
colkey+         length = 0.3, width = 0.5, col.axis = "white", 
colkey+         col.ticks = "white", at  = c(100, 140, 180), 
colkey+         labels = c("a", "b", "c"), font = 2)

colkey>  colkey(clim = range(volcano), dist = -0.3, shift = -0.2,
colkey+         side = 2, add = TRUE, clab = "key 4", col.clab = "white",
colkey+         length = 0.3, width = 0.5, col.axis = "white", 
colkey+         col.ticks = "white", col.box = "red", cex.axis = 0.8,
colkey+         las = 3)

colkey> ## =======================================================================
colkey> ##  colorkey in other plots
colkey> ## =======================================================================
colkey> 
colkey>  par(mfrow = c(1, 1))

colkey>  par(mar = par("mar") + c(0, 0, -2, 0))

colkey>  image2D(volcano, clab = "height, m", 
colkey+        colkey = list(dist = -0.15, shift = 0.2,
colkey+        side = 3, length = 0.5, width = 0.5, line.clab = 2.5,
colkey+        cex.clab = 2, col.clab = "white", col.axis = "white", 
colkey+        col.ticks = "white", cex.axis = 0.8))

colkey> ## =======================================================================
colkey> ## Several color keys in composite plot
colkey> ## =======================================================================
colkey> 
colkey>  persp3D(z = volcano, zlim = c(-60, 200), phi = 20, bty = "b",   
colkey+     colkey = list(length = 0.2, width = 0.4, shift = 0.15,
colkey+       cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90,
colkey+     clab = c("height","m"), plot = FALSE)

colkey> # create gradient in x-direction
colkey>  Vx <- volcano[-1, ] - volcano[-nrow(volcano), ]

colkey> # add as image with own color key, at bottom 
colkey>  image3D(z = -60, colvar = Vx/10, add = TRUE, 
colkey+     colkey = list(length = 0.2, width = 0.4, shift = -0.15,
colkey+       cex.axis = 0.8, cex.clab = 0.85),
colkey+    clab = c("gradient","m/m"), plot = TRUE)

colkey> ## =======================================================================
colkey> ## categorical colors; use addlines = TRUE to separate colors
colkey> ## =======================================================================
colkey> 
colkey>  with(iris, scatter3D(x = Sepal.Length, y = Sepal.Width, 
colkey+    z = Petal.Length, colvar = as.integer(Species), 
colkey+    col = c("orange", "green", "lightblue"), pch = 16, cex = 2,  
colkey+    clim = c(1, 3), ticktype = "detailed", phi = 20,
colkey+    xlab = "Sepal Length", ylab = "Sepal Width", 
colkey+    zlab = "Petal Length",  main = "iris",
colkey+    colkey = list(at = c(1.33, 2, 2.66), side = 1, 
colkey+    addlines = TRUE, length = 0.5, width = 0.5,
colkey+    labels = c("setosa", "versicolor", "virginica") )))

colkey> # reset plotting parameters
colkey>  par(mfrow = pm)

colkey>  par(mar = pmar)

jet.cl> # save plotting parameters
jet.cl>  pm <- par("mfrow")

jet.cl>  pmar <- par("mar")

jet.cl> ## =======================================================================
jet.cl> ## Transparency and various color schemes
jet.cl> ## =======================================================================
jet.cl> 
jet.cl>  par(mfrow = c(3, 3))

jet.cl>  for (alph in c(0.25, 0.75))
jet.cl+    image2D(volcano, alpha = alph, 
jet.cl+          main = paste("jet.col, alpha = ", alph))  

jet.cl>  image2D(volcano, main = "jet.col")

jet.cl>  image2D(volcano, col = jet2.col(100), main = "jet2.col")

jet.cl>  image2D(volcano, col = gg.col(100), main = "gg.col")

jet.cl>  image2D(volcano, col = gg2.col(100), main = "gg2.col")

jet.cl>  image2D(volcano, col = rainbow(100), main = "rainbow")

jet.cl>  image2D(volcano, col = terrain.colors(100), main = "terrain.colors")

jet.cl>  image2D(volcano, col = ramp.col(c("blue", "yellow", "green", "red")),
jet.cl+        main = "ramp.col")  

jet.cl> ## =======================================================================
jet.cl> ## Shading, lighting -  one color
jet.cl> ## =======================================================================
jet.cl> 
jet.cl> # create grid matrices
jet.cl>  X      <- seq(0, pi, length.out = 50)

jet.cl>  Y      <- seq(0, 2*pi, length.out = 50)

jet.cl>  M      <- mesh(X, Y)

jet.cl>  phi    <- M$x

jet.cl>  theta  <- M$y

jet.cl> # x, y and z grids
jet.cl>  x <- sin(phi) * cos(theta)

jet.cl>  y <- cos(phi)

jet.cl>  z <- sin(phi) * sin(theta)

jet.cl> # these are the defaults
jet.cl>  p <- list(ambient = 0.3, diffuse = 0.6, specular = 1.,
jet.cl+            exponent = 20, sr = 0, alpha = 1)

jet.cl>  par(mfrow = c(3, 3), mar = c(0, 0, 0, 0))

jet.cl>  Col <- "red"

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, shade = 0.9)

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = TRUE)  

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(ambient = 0))

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(diffuse = 0))

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(diffuse = 1))

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(specular = 0))

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(exponent = 5))

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(exponent = 50))

jet.cl>  surf3D(x, y, z, box = FALSE, col = Col, lighting = list(sr = 1)) 

jet.cl> ## =======================================================================
jet.cl> ## Shading, lighting with default colors
jet.cl> ## =======================================================================
jet.cl> 
jet.cl>  x <- seq(-pi, pi, len = 100)

jet.cl>  y <- seq(-pi, pi, len = 100)

jet.cl>  grid <- mesh(x, y)

jet.cl>  z    <- with(grid, cos(x) * sin(y))

jet.cl>  cv   <- with(grid, -cos(y) * sin(x))

jet.cl> # lphi = 180, ltheta = -130  - good for shade
jet.cl> # lphi = 90, ltheta = 0  - good for lighting
jet.cl> 
jet.cl>  par(mfrow = c(2, 2))

jet.cl>  persp3D(z = z, x = x, y = y, colvar = cv, zlim = c(-3, 3), colkey = FALSE)

jet.cl>  persp3D(z = z, x = x, y = y, colvar = cv, zlim = c(-3, 3), 
jet.cl+        lighting = TRUE, colkey = FALSE)

jet.cl>  persp3D(z = z, x = x, y = y, colvar = cv, zlim = c(-3, 3), 
jet.cl+        shade = 0.25, colkey = FALSE)

jet.cl>  persp3D(z = z, x = x, y = y, colvar = cv, zlim = c(-3, 3), 
jet.cl+        lighting = TRUE, lphi = 90, ltheta = 0, colkey = FALSE)

jet.cl> ## =======================================================================
jet.cl> ## transparency of a vector of colors
jet.cl> ## =======================================================================
jet.cl> 
jet.cl>  par(mfrow = c(1, 1))

jet.cl>  x <- runif(19)

jet.cl>  y <- runif(19)

jet.cl>  z <- runif(19)

jet.cl> # split into 5 sections (polygons) 
jet.cl>  ii <- seq(4, 19, by = 4)

jet.cl>  x[ii] <- y[ii] <- z[ii] <- NA

jet.cl>  polygon3D(x, y, z, border = "black", lwd = 2, 
jet.cl+    col = alpha.col(c("red", "lightblue", "yellow", "green", "black"), 
jet.cl+                   alpha = 0.4))

jet.cl> # the same, now passing alpha as an argument to polygon3D:
jet.cl> ## Not run: 
jet.cl> ##D  polygon3D(x, y, z, border = "black", lwd = 2, 
jet.cl> ##D    col = c("red", "lightblue", "yellow", "green", "black"), 
jet.cl> ##D                   alpha = 0.4)
jet.cl> ## End(Not run)
jet.cl> # reset plotting parameters
jet.cl>  par(mfrow = pm)

jet.cl>  par(mar = pmar)

prspbx> # save plotting parameters                            
prspbx>  pm   <- par("mfrow")

prspbx>  pmar <- par("mar")

prspbx> ## ========================================================================
prspbx> ## The 4 predefined box types
prspbx> ## ========================================================================
prspbx> 
prspbx>  par(mfrow = c(2, 2), mar = c(1, 1, 1, 1))

prspbx> # box type with only backward panels
prspbx>  perspbox(z = volcano, bty = "b", ticktype = "detailed", d = 2, 
prspbx+           main  = "bty = 'b'")

prspbx> # box as in 'persp'
prspbx>  perspbox(z = volcano, bty = "f", ticktype = "detailed", 
prspbx+           d = 2, main  = "bty = 'f'")

prspbx> # back panels with gridlines, detailed axes
prspbx>  perspbox(z = volcano, bty = "b2", ticktype = "detailed", 
prspbx+           d = 2, main  = "bty = 'b2'")

prspbx> # ggplot-type, simple axes 
prspbx>  perspbox(z = volcano, bty = "g", 
prspbx+           d = 2, main  = "bty = 'g'")

prspbx> ## ========================================================================
prspbx> ## A user-defined box
prspbx> ## ========================================================================
prspbx> 
prspbx>  par(mfrow = c(1, 1))

prspbx>  perspbox(z = diag(2), bty = "u", ticktype = "detailed", 
prspbx+           col.panel = "gold", col.axis = "white",  
prspbx+           scale = FALSE, expand = 0.4, 
prspbx+           col.grid = "grey", main = "user-defined")

prspbx> # restore plotting parameters
prspbx>  par(mfrow = pm)

prspbx>  par(mar = pmar)

mesh> ## ========================================================================
mesh> ## 2-D mesh
mesh> ## ========================================================================
mesh> 
mesh>  x <- c(-1 , 0, 1)

mesh>  y <- 1 : 4

mesh> # 2-D mesh
mesh>  (M <- mesh(x, y))
$x
     [,1] [,2] [,3] [,4]
[1,]   -1   -1   -1   -1
[2,]    0    0    0    0
[3,]    1    1    1    1

$y
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4


mesh> # calculate with this mesh
mesh>  V <- with (M, x/2 * sin(y))

mesh> # same as:
mesh>  V2 <- outer(x, y, FUN = function(x, y) x/2*sin(y))

mesh> ## ========================================================================
mesh> ## 3-D mesh
mesh> ## ========================================================================
mesh> 
mesh>  x <- y <- z <- c(-1 , 0, 1)

mesh> # 3-D mesh
mesh>  (M <- mesh(x, y, z))
$x
, , 1

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1

, , 2

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1

, , 3

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1


$y
, , 1

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1

, , 2

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1

, , 3

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1


$z
, , 1

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]   -1   -1   -1
[3,]   -1   -1   -1

, , 2

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

, , 3

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1



mesh> # calculate with 3-D mesh
mesh>  V <- with (M, x/2 * sin(y) *sqrt(z+2))

mesh> # plot result
mesh>  scatter3D(M$x, M$y, M$z, V, pch = "+", cex = 3, colkey = FALSE)

trns3D> ## ========================================================================
trns3D> ## 3-D mesh
trns3D> ## ========================================================================
trns3D> 
trns3D>  x <- y <- z <- c(-1 , 0, 1)

trns3D> # plot a 3-D mesh
trns3D>  (M <- mesh(x, y, z))
$x
, , 1

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1

, , 2

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1

, , 3

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    0    0    0
[3,]    1    1    1


$y
, , 1

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1

, , 2

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1

, , 3

     [,1] [,2] [,3]
[1,]   -1    0    1
[2,]   -1    0    1
[3,]   -1    0    1


$z
, , 1

     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]   -1   -1   -1
[3,]   -1   -1   -1

, , 2

     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0

, , 3

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1



trns3D> # plot result
trns3D>  pmat <- scatter3D(M$x, M$y, M$z, pch = "+", cex = 3, colkey = FALSE)

trns3D> # add line
trns3D>  XY <- trans3D(x = c(-1, 1), y = c(-1, 1), z = c(-1, 1), pmat = pmat) 

trns3D>  lines(XY, lwd = 2, col = "blue")

trns3D> ## ========================================================================
trns3D> ## Example 2
trns3D> ## ========================================================================
trns3D> 
trns3D>  pmat <- perspbox (z = diag(2))

trns3D>  XY <- trans3D(x = runif(30), y = runif(30), z = runif(30), pmat = pmat) 

trns3D>  polygon(XY, col = "darkblue")

plt.pl> # save plotting parameters                            
plt.pl>  pm   <- par("mfrow")

plt.pl>  pmar <- par("mar")

plt.pl> ## ========================================================================
plt.pl> ## The volcano
plt.pl> ## ========================================================================
plt.pl> 
plt.pl>  par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))

plt.pl> # The volcano at lower resolution
plt.pl>  x <- seq(1, nrow(volcano), by = 2)

plt.pl>  y <- seq(1, ncol(volcano), by = 2)

plt.pl>  V <- volcano[x,y]

plt.pl>  persp3D(z = V)

plt.pl> # rotate
plt.pl>  plotdev(theta = 0)

plt.pl> # light and transparence
plt.pl>  plotdev(lighting  = TRUE, lphi = 90, alpha = 0.6)  

plt.pl> # zoom
plt.pl>  plotdev(xlim = c(0.2, 0.6), ylim = c(0.2, 0.6), phi = 60) 

plt.pl> ## ========================================================================
plt.pl> ## Two spheres 
plt.pl> ## ========================================================================
plt.pl> 
plt.pl>  par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))

plt.pl> # create a sphere
plt.pl>  M  <- mesh(seq(0, 2*pi, length.out = 30),
plt.pl+             seq(0,   pi, length.out = 30))

plt.pl>  u  <- M$x ; v  <- M$y

plt.pl>  x <- cos(u)*sin(v)

plt.pl>  y <- sin(u)*sin(v)

plt.pl>  z <- cos(v)

plt.pl>  surf3D(x = 2*x, y = 2*y, z = 2*z, 
plt.pl+         colvar = NULL, lighting = TRUE, #plot = FALSE,
plt.pl+         facets = NA, col = "blue", lwd = 5)

plt.pl>  surf3D(x, y, z, colvar = NULL, lighting = TRUE, 
plt.pl+         col = "red", add = TRUE)

plt.pl>  names(getplist())
 [1] "type"     "plt"      "persp"    "alpha"    "setlim"   "xlim"    
 [7] "ylim"     "zlim"     "scalefac" "mat"      "dot"      "imgnr"   
[13] "img"      "poly"    

plt.pl> # plot with different view:
plt.pl>  plotdev(phi = 0)  

plt.pl> ## Not run: 
plt.pl> ##D   # will plot same 3-D graph to pdf
plt.pl> ##D  pdf(file = "save.pdf")
plt.pl> ##D  plotdev()
plt.pl> ##D  dev.off()
plt.pl> ## End(Not run)
plt.pl>              
plt.pl> ## ========================================================================
plt.pl> ## Two spheres and two planes 
plt.pl> ## ========================================================================
plt.pl> 
plt.pl>  par(mar = c(2, 2, 2, 2))

plt.pl> # equation of a sphere
plt.pl>  M  <- mesh(seq(0, 2*pi, length.out = 100),                                     -
plt.pl+             seq(0,   pi, length.out = 100))

plt.pl>  u  <- M$x ; v  <- M$y

plt.pl>  x <- cos(u)*sin(v)

plt.pl>  y <- sin(u)*sin(v)

plt.pl>  z <- cos(v)

plt.pl>  surf3D(x, y, z, colvar = z, 
plt.pl+         theta = 45, phi = 20, bty = "b",
plt.pl+         xlim = c(-1.5, 1.5), ylim = c(-1, 2), 
plt.pl+         zlim = c(-1.5, 1.5), plot = FALSE)

plt.pl> # add a second sphere, shifted 1 unit to the right on y-axis; 
plt.pl> # no facets drawn for this sphere 
plt.pl>  surf3D (x, y+1, z, colvar = z, add = TRUE, 
plt.pl+          facets = FALSE, plot = FALSE)

plt.pl> # define a plane at z = 0
plt.pl>  Nx <- 100

plt.pl>  Ny <- 100

plt.pl>  x <- seq(-1.5, 1.5, length.out = Nx)

plt.pl>  y <- seq(-1, 2, length.out = Ny)

plt.pl>  image3D (x = x, y = y, z = 0, add = TRUE, colvar = NULL, 
plt.pl+           col = "blue", facets = TRUE, plot = FALSE)

plt.pl> # another, small plane at y = 0 - here x and y have to be matrices!
plt.pl>  x <- seq(-1., 1., length.out = 50)

plt.pl>  z <- seq(-1., 1., length.out = 50)

plt.pl>  image3D (x = x, y = 0, z = z, colvar = NULL, 
plt.pl+          add = TRUE, col = NA, border = "blue", 
plt.pl+          facets = TRUE, plot = TRUE)       

plt.pl> ## Not run: 
plt.pl> ##D   # rotate 
plt.pl> ##D  for (angle in seq(0, 360, by = 10)) 
plt.pl> ##D    plotdev(theta = angle)
plt.pl> ## End(Not run)
plt.pl> 
plt.pl> ## ========================================================================
plt.pl> ## Zooming, rescaling, lighting,...
plt.pl> ## ========================================================================
plt.pl> 
plt.pl>  par(mfrow = c(2, 2)) 

plt.pl> # The volcano
plt.pl>  x <- seq(1, nrow(volcano), by = 2)

plt.pl>  y <- seq(1, ncol(volcano), by = 2)

plt.pl>  V <- volcano[x,y]

plt.pl> # plot the volcano
plt.pl>  persp3D (x, y, z = V, colvar = V, theta = 10, phi = 20, 
plt.pl+           box = FALSE, scale = FALSE, expand = 0.3, 
plt.pl+           clim = range(V), plot = FALSE)

plt.pl> # add a plane (image) at z = 170; jetcolored, transparant: only border
plt.pl>  image3D(x, y, z = 170, add = TRUE, clim = range(V), 
plt.pl+          colvar = V, facets = NA, plot = FALSE, colkey = FALSE)

plt.pl> # add a contour (image) at z = 170; jetcolored, 
plt.pl>  contour3D(x, y, z = 170, add = TRUE, clim = range(V),
plt.pl+            colvar = V, plot = FALSE, colkey = FALSE)

plt.pl> # plot it  - 
plt.pl>  plot(getplist())   #  same as plotdev()

plt.pl> # plot but with different expansion
plt.pl>  plotdev(expand = 1)

plt.pl> # other perspective, and shading
plt.pl>  plotdev(d = 2, r = 10, shade = 0.3)

plt.pl> # zoom and rotate
plt.pl>  plotdev(xlim = c(10, 30), ylim = c(20, 30), phi = 50)

plt.pl> ## ========================================================================
plt.pl> ## Using setplist
plt.pl> ## ========================================================================
plt.pl> 
plt.pl>  polygon3D(runif(3), runif(3), runif(3))

plt.pl> # retrieve plotting list
plt.pl>  plist <- getplist()

plt.pl>  names(plist)
 [1] "type"     "plt"      "persp"    "alpha"    "setlim"   "xlim"    
 [7] "ylim"     "zlim"     "scalefac" "mat"      "dot"      "imgnr"   
[13] "img"      "poly"    

plt.pl>  plist$poly
$x
          [,1]
[1,] 0.3147765
[2,] 0.5926986
[3,] 0.6017956
[4,]        NA

$y
           [,1]
[1,] 0.60047536
[2,] 0.82227459
[3,] 0.04894024
[4,]         NA

$z
          [,1]
[1,] 0.6877397
[2,] 0.7195563
[3,] 0.8633404
[4,]        NA

$col
[1] "grey"

$border
[1] NA

$lwd
[1] 1

$lty
[1] 1

$alpha
[1] NA

$isimg
[1] 0

$proj
[1] -9.339463

attr(,"class")
[1] "poly"

plt.pl> # change copy of plotting list
plt.pl>  plist$poly$col <- "red"

plt.pl> # update internal plotting list
plt.pl>  setplist(plist)

plt.pl> # plot updated list
plt.pl>  plotdev()

plt.pl> ## ========================================================================
plt.pl> ## Using selectplist
plt.pl> ## ========================================================================
plt.pl> 
plt.pl>  polygon3D(runif(10), runif(10), runif(10), col = "red", 
plt.pl+    alpha = 0.2, plot = FALSE, ticktype = "detailed", 
plt.pl+    xlim = c(0,1), ylim = c(0, 1), zlim = c(0, 1))

plt.pl>  polygon3D(runif(10)*0.5, runif(10), runif(10), col = "yellow", 
plt.pl+    alpha = 0.2, plot = FALSE, add = TRUE)

plt.pl>  polygon3D(runif(10)*0.5+0.5, runif(10), runif(10), col = "green", 
plt.pl+    alpha = 0.2, plot = FALSE, add = TRUE)

plt.pl>  points3D(runif(10), runif(10), runif(10), col = "blue", 
plt.pl+    add = TRUE, plot = FALSE)

plt.pl>  segments3D(x0 = runif(10), y0 = runif(10), z0 = runif(10), 
plt.pl+    x1 = runif(10), y1 = runif(10), z1 = runif(10), 
plt.pl+    colvar = 1:10, add = TRUE, lwd = 3)

plt.pl> # retrieve plotting list
plt.pl>  plist <- getplist()

plt.pl> # selection function 
plt.pl>  SS <- function (x, y, z)  {
plt.pl+    sel <- rep(TRUE, length.out = length(x))
plt.pl+    sel[x < 0.5] <- FALSE
plt.pl+    return(sel)
plt.pl+  } 

plt.pl> # The whole polygon will be removed or kept.  
plt.pl>  plot(x = selectplist(plist, SS), 
plt.pl+    xlim = c(0, 1), ylim = c(0, 1), zlim = c(0, 1))

plt.pl> # restore plotting parameters
plt.pl>  par(mfrow = pm)

plt.pl>  par(mar = pmar)

ImgOcn> # save plotting parameters
ImgOcn>  pm <- par("mfrow")

ImgOcn>  mar <- par("mar")

ImgOcn> ## =======================================================================
ImgOcn> ## Images of the hypsometry
ImgOcn> ## =======================================================================
ImgOcn> 
ImgOcn>  par(mfrow = c(2, 2))

ImgOcn>  image2D(Hypsometry, asp = TRUE, xlab = expression(degree*E), 
ImgOcn+    ylab =  expression(degree*N), contour = TRUE)

ImgOcn> # remove ocean
ImgOcn>  zz         <- Hypsometry$z

ImgOcn>  zz[zz < 0] <- NA

ImgOcn>  image2D(zz, x = Hypsometry$x, y = Hypsometry$y, NAcol = "black")

ImgOcn> ## =======================================================================
ImgOcn> ## A short version for plotting the Ocean's bathymetry
ImgOcn> ## =======================================================================
ImgOcn> 
ImgOcn>  ImageOcean(asp = TRUE, contour = TRUE)

ImgOcn>  ImageOcean(col = "white", 
ImgOcn+    contour = list(levels = seq(-6000, 0, by = 2000)))

ImgOcn> ## =======================================================================
ImgOcn> ## A complex image of part of the ocean 
ImgOcn> ## =======================================================================
ImgOcn> 
ImgOcn> # elaborate version
ImgOcn>  par(mfrow = c(1, 1), mar = c(2, 2, 2, 2))

ImgOcn>  ii <- which(Hypsometry$x > -50 & Hypsometry$x < -20)

ImgOcn>  jj <- which(Hypsometry$y >  10 & Hypsometry$y <  40)

ImgOcn> # Draw empty persp box
ImgOcn>  zlim <- c(-10000, 0)

ImgOcn>  pmat <- perspbox(z = Hypsometry$z[ii, jj], 
ImgOcn+                   xlab = "longitude", ylab = "latitude", zlab = "depth", 
ImgOcn+                   expand = 0.5, d = 2, zlim = zlim, phi = 20, theta = 30,
ImgOcn+                   colkey = list(side = 1))

ImgOcn> # A function that makes a black panel with grey edge:
ImgOcn>  panelfunc <- function(x, y, z) {
ImgOcn+     XY <- trans3D(x, y, z, pmat = pmat)
ImgOcn+     polygon(XY$x, XY$y, col = "black", border = "grey")
ImgOcn+  }

ImgOcn> # left panel
ImgOcn>  panelfunc(x = c(0, 0, 0, 0), y = c(0, 0, 1, 1), 
ImgOcn+            z = c(zlim[1], zlim[2], zlim[2], zlim[1]))

ImgOcn> # back panel
ImgOcn>  panelfunc(x = c(0, 0, 1, 1), y = c(1, 1, 1, 1),
ImgOcn+            z = c(zlim[1], zlim[2], zlim[2], zlim[1]))

ImgOcn> # bottom panel
ImgOcn>  panelfunc(x = c(0, 0, 1, 1), y = c(0, 1, 1, 0),
ImgOcn+            z = c(zlim[1], zlim[1], zlim[1], zlim[1]))

ImgOcn> # Actual bathymetry, 2 times increased resolution, with contours
ImgOcn>  persp3D(z = Hypsometry$z[ii,jj], add = TRUE, resfac = 2,  
ImgOcn+        contour = list(col = "grey", side = c("zmin", "z")), 
ImgOcn+        zlim = zlim, clab = "depth, m", 
ImgOcn+        colkey = list(side = 1, length = 0.5, dist = -0.1))

ImgOcn> # shorter alternative for same plot, higher resolution
ImgOcn> ## Not run: 
ImgOcn> ##D  persp3D(z = Hypsometry$z[ii,jj], resfac = 4,  
ImgOcn> ##D        contour = list(col = "grey", side = c("zmin", "z")), 
ImgOcn> ##D        zlim = zlim, clab = "depth, m", bty = "bl2",
ImgOcn> ##D        xlab = "longitude", ylab = "latitude", zlab = "depth", 
ImgOcn> ##D        expand = 0.5, d = 2, phi = 20, theta = 30,
ImgOcn> ##D        colkey = list(side = 1, length = 0.5, dist = -0.1))
ImgOcn> ## End(Not run)
ImgOcn> 
ImgOcn> # reset plotting parameters
ImgOcn>  par(mfrow = pm)

ImgOcn>  par(mar = mar)

Oxsat> # save plotting parameters
Oxsat>  pm <- par("mfrow")

Oxsat> ## ========================================================================
Oxsat> ## plot all surface data
Oxsat> ## ========================================================================
Oxsat> 
Oxsat>  par(mfrow = c(1, 1))

Oxsat>  image2D(z = Oxsat$val[ , , 1], x = Oxsat$lon, y = Oxsat$lat,
Oxsat+        main = "surface oxygen saturation (%) for 2005")

Oxsat> ## ========================================================================
Oxsat> ## plot a selection of latitude-depth profiles; input is an array
Oxsat> ## ========================================================================
Oxsat> 
Oxsat>  lon <- Oxsat$lon

Oxsat>  image2D (z = Oxsat$val, margin = c(2, 3), x = Oxsat$lat, 
Oxsat+         y = Oxsat$depth, subset = (lon > 18 & lon < 23),
Oxsat+         ylim = c(5500, 0), NAcol = "black", zlim = c(0, 110),
Oxsat+         xlab = "latitude", ylab = "depth, m")

Oxsat>  ImageOcean()

Oxsat>  abline ( v = lon[lon > 18 & lon < 23])

Oxsat> ## ========================================================================
Oxsat> ## plot with slices
Oxsat> ## ========================================================================
Oxsat> 
Oxsat>  par(mfrow = c(1, 1))

Oxsat>  ii <- which (Oxsat$lon > -90 & Oxsat$lon < 90)

Oxsat>  jj <- which (Oxsat$lat > 0 & Oxsat$lat < 90)

Oxsat>  xs <- Oxsat$lon[ii[length(ii)]]  # E boundary

Oxsat>  ys <- Oxsat$lat[jj[1]]           # S boundary

Oxsat>  slice3D(colvar = Oxsat$val[ii,jj,], x = Oxsat$lon[ii],  
Oxsat+         y = Oxsat$lat[jj], z = -Oxsat$depth,
Oxsat+         NAcol = "black", xs = xs, ys = ys, zs = 0, 
Oxsat+         theta = 35, phi = 50, colkey = list(length = 0.5),
Oxsat+         expand = 0.5, ticktype = "detailed",
Oxsat+         clab = "%", main = "Oxygen saturation", 
Oxsat+         xlab = "longitude", ylab = "latitude", zlab = "depth")

Oxsat> # restore plotting parameters
Oxsat>  par(mfrow = pm)

plot3D documentation built on May 22, 2021, 5:06 p.m.

Related to plot3D in plot3D...