-
-
Notifications
You must be signed in to change notification settings - Fork 5
/
persp.R
114 lines (77 loc) · 3.13 KB
/
persp.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
> ### Demos for persp() plots -- things not in example(persp)
> ### -------------------------
>
> require(datasets)
> require(grDevices); require(graphics)
> ## (1) The Obligatory Mathematical surface.
> ## Rotated sinc function.
>
> x <- seq(-10, 10, length.out = 50)
> y <- x
> rotsinc <- function(x,y)
+ {
+ sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }
+ 10 * sinc( sqrt(x^2+y^2) )
+ }
> sinc.exp <- expression(z == Sinc(sqrt(x^2 + y^2)))
> z <- outer(x, y, rotsinc)
> oldpar <- par(bg = "white")
> persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")
Hit <Return> to see next plot:
> title(sub=".")## work around persp+plotmath bug
> title(main = sinc.exp)
> persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
+ ltheta = 120, shade = 0.75, ticktype = "detailed",
+ xlab = "X", ylab = "Y", zlab = "Z")
Hit <Return> to see next plot:
> title(sub=".")## work around persp+plotmath bug
> title(main = sinc.exp)
> ## (2) Visualizing a simple DEM model
>
> z <- 2 * volcano # Exaggerate the relief
> x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)
> y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)
> persp(x, y, z, theta = 120, phi = 15, scale = FALSE, axes = FALSE)
Hit <Return> to see next plot:
> ## (3) Now something more complex
> ## We border the surface, to make it more "slice like"
> ## and color the top and sides of the surface differently.
>
> z0 <- min(z) - 20
> z <- rbind(z0, cbind(z0, z, z0), z0)
> x <- c(min(x) - 1e-10, x, max(x) + 1e-10)
> y <- c(min(y) - 1e-10, y, max(y) + 1e-10)
> fill <- matrix("green3", nrow = nrow(z)-1, ncol = ncol(z)-1)
> fill[ , i2 <- c(1,ncol(fill))] <- "gray"
> fill[i1 <- c(1,nrow(fill)) , ] <- "gray"
> par(bg = "lightblue")
> persp(x, y, z, theta = 120, phi = 15, col = fill, scale = FALSE, axes = FALSE)
Hit <Return> to see next plot:
> title(main = "Maunga Whau\nOne of 50 Volcanoes in the Auckland Region.",
+ font.main = 4)
> par(bg = "slategray")
> persp(x, y, z, theta = 135, phi = 30, col = fill, scale = FALSE,
+ ltheta = -120, lphi = 15, shade = 0.65, axes = FALSE)
Hit <Return> to see next plot:
> ## Don't draw the grid lines : border = NA
> persp(x, y, z, theta = 135, phi = 30, col = "green3", scale = FALSE,
+ ltheta = -120, shade = 0.75, border = NA, box = FALSE)
Hit <Return> to see next plot:
> ## `color gradient in the soil' :
> fcol <- fill ; fcol[] <- terrain.colors(nrow(fcol))
> persp(x, y, z, theta = 135, phi = 30, col = fcol, scale = FALSE,
+ ltheta = -120, shade = 0.3, border = NA, box = FALSE)
Hit <Return> to see next plot:
> ## `image like' colors on top :
> fcol <- fill
> zi <- volcano[ -1,-1] + volcano[ -1,-61] +
+ volcano[-87,-1] + volcano[-87,-61] ## / 4
> fcol[-i1,-i2] <-
+ terrain.colors(20)[cut(zi,
+ stats::quantile(zi, seq(0,1, length.out = 21)),
+ include.lowest = TRUE)]
> persp(x, y, 2*z, theta = 110, phi = 40, col = fcol, scale = FALSE,
+ ltheta = -120, shade = 0.4, border = NA, box = FALSE)
Hit <Return> to see next plot:
> ## reset par():
> par(oldpar)