mark 学习
params <- function(N, breaks, p=seq(0.001, 1, length=100)) {
list(N=N, T=1/breaks, p=p, q=1-p)
}
pdfcomp <- function(comp, params) {
n <- params$T
p <- params$p
q <- params$q
y <- round(comp/n)
choose(n, comp)*p^comp*q^(n-comp) / (1 - q^n)
}
# Expected num sherds (for a vessel) [=completeness]
expcomp <- function(params) {
params$T*params$p/(1-params$q^params$T)
}
# Variance of num sherds (for a vessel)
varcomp <- function(params) {
n <- params$T
p <- params$p
q <- params$q
# From Johnson & Kotz
(n*p*q / (1 - q^n)) - (n^2*p^2*q^n / (1 - q^n)^2)
# n^2 times Thomas Yee's formula
# n^2*((p*(1 + p*(n - 1)) / (n*(1 - q^n))) - (p^2 / (1 - q^n)^2))
}
# Expected value of completeness (for a sample of vessels)
expmeancomp <- function(params) {
expcomp(params)
}
# Variance of completeness (for a sample of vessels)
# Use the expected number of vessels in sample as denominator
varmeancomp <- function(params) {
varcomp(params)/(numvess(params))
}
numvess <- function(params) {
params$N*(1-params$q^params$T)
}
ecomp <- function(p, T, comp) {
q <- 1 - p
T*p/(1 - q^T) - comp
}
estN <- function(comp, broke, n) {
T <- 1/broke
n / (1 - (1 - uniroot(ecomp, c(0.00001, 1), T=T, comp=comp)$root)^T)
}
nvessscale <- function(params, xlim, ylim, new=TRUE) {
if (new)
par(new=TRUE)
plot(0:1, c(1, params$N), type="n", axes=!new, ann=FALSE,
xlim=xlim, ylim=ylim)
}
compscale <- function(params, xlim, ylim, new=TRUE) {
if (new)
par(new=TRUE)
plot(0:1, c(1, params$T), type="n", axes=!new, ann=FALSE,
xlim=xlim, ylim=ylim)
}
lowerCI <- function(p, N, breaks, lb) {
params <- params(N, breaks, p)
expmeancomp(params) - 2*sqrt(varmeancomp(params)) - lb
}
upperCI <- function(p, N, breaks, lb) {
params <- params(N, breaks, p)
expmeancomp(params) + 2*sqrt(varmeancomp(params)) - lb
}
critP <- function(comp, params) {
c(uniroot(lowerCI, c(0.00001, 1), N=params$N,
breaks=1/params$T, lb=max(comp))$root,
if (upperCI(0.00001, params$N, 1/params$T, min(comp)) > 0) 0
else uniroot(upperCI, c(0.00001, 1), N=params$N,
breaks=1/params$T, lb=min(comp))$root)
}
anncomp <- function(params, comp, xlim, ylim, cylim) {
cp <- critP(comp, params)
nv <- numvess(params(params$N, 1/params$T, cp))
nvessscale(params, xlim, ylim)
polygon(c(cp[2], cp[2], 0, 0, cp[1], cp[1]),
c(0, nv[2], nv[2], nv[1], nv[1], 0),
col="grey90", border=NA)
text(0, nv[1], paste(round(nv[1]),
" (", round(100*nv[1]/params$N), "%)", sep=""),
adj=c(0, 0), col="grey")
text(0, nv[2], paste(round(nv[2]),
" (", round(100*nv[2]/params$N), "%)", sep=""),
adj=c(0, 1), col="grey")
compscale(params, xlim, cylim)
segments(1, min(comp), cp[2], comp, col="grey")
segments(1, max(comp), cp[1], comp, col="grey")
text(1, comp, paste(comp, collapse="-"), adj=c(1, 0), col="grey")
}
plotPars <- function(params, comp, xlim=NULL, ylim=NULL) {
mean <- expmeancomp(params)
var <- 2*sqrt(varmeancomp(params))
lb <- mean - var
ub <- mean + var
par(mar=c(5, 4, 4, 4))
if (is.null(ylim))
cylim <- ylim
else
cylim <- c(1 + ((ylim[1] - 1)/(params$N - 1))*(params$T - 1),
1 + ((ylim[2] - 1)/(params$N - 1))*(params$T - 1))
nvessscale(params, xlim, ylim, new=FALSE)
compscale(params, xlim, cylim)
polygon(c(params$p, rev(params$p)), c(lb, rev(ub)),
col="grey90", border=NA)
anncomp(params, comp, xlim, ylim, cylim)
nvessscale(params, xlim, ylim)
mtext("Number of Vessels", side=2, line=3)
mtext("Sampling Fraction", side=1, line=3)
lines(params$p, numvess(params))
par(new=TRUE)
compscale(params, xlim, cylim)
mtext("Completeness", side=4, line=3)
axis(4)
lines(params$p, mean, lty="dashed")
lines(params$p, lb, lty="dotted")
lines(params$p, ub, lty="dotted")
mtext(paste("N = ", round(params$N),
" brokenness = ", round(1/params$T, 3), sep=""),
side=3, line=2)
}
par(cex=0.8, mar=c(3, 3, 3, 3),bg=7)
p6 <- params(estN(1.2, 0.5, 200), 0.5)
plotPars(p6, 1.2)
nvessscale(p6, NULL, NULL)
pcrit <- 1 - (1 - 200/estN(1.2, 0.5, 200))^(1/p6$T)
lines(c(0, pcrit), c(200, 200))
lines(c(pcrit, pcrit), c(200, 0))
##############################################################################
#
# Comment:
#
# An example of a one-off image drawn using the grid system.
#
# The code is somewhat modular and general, with functions
# for producing different shapes, but the sizes and
# locations used in this particular image assume a 2:1 aspect ratio.
#
# The gradient-fill background (dark at the top to lighter at the
# bottom) is achieved by filling multiple overlapping polygons with
# slowly changing shades of grey.
#
install.packages("RGraphics")
library(RGraphics)
vignette("RGraphics",package="RGraphics")
pushViewport(viewport(xscale=c(0, 1), yscale=c(0.5, 1),
clip=TRUE))
res <- 50
for (i in 1:res)
grid.rect(y=1 - (i-1)/res, just="top",
gp=gpar(col=NULL, fill=grey(0.5*i/res)))
moon <- function(x, y, size) {
angle <- seq(-90, 90, length=50)/180*pi
x1 <- x + size*cos(angle)
y1 <- y + size*sin(angle)
mod <- 0.7
x2 <- x + mod*(x1 - x)
grid.polygon(c(x1, rev(x2)), c(y1, rev(y1)),
default.unit="native",
gp=gpar(col=7, fill="white"))
}
moon(.1, .9, .02)
star <- function(x, y, size) {
x1 <- c(x, x + size*.1, x + size*.5, x + size*.1,
x, x - size*.1, x - size*.5, x - size*.1) + .05
y1 <- c(y - size, y - size*.1, y, y + size*.1,
y + size*.7, y + size*.1, y, y - size*.1) + .05
grid.polygon(x1, y1,
default.unit="native",
gp=gpar(col=5, fill="white"))
}
star(.5, .7, .02)
star(.8, .9, .02)
star(.72, .74, .02)
star(.62, .88, .02)
grid.circle(runif(20, .2, 1), runif(20, .6, 1), r=.002,
default.unit="native",
gp=gpar(col=5, fill="white"))
hill <- function(height=0.1, col=4) {
n <- 100
x <- seq(0, 1, length=n)
y1 <- sin(runif(1) + x*2*pi)
y2 <- sin(runif(1) + x*4*pi)
y3 <- sin(runif(1) + x*8*pi)
y <- 0.6 + height*((y1 + y2 + y3)/3)
grid.polygon(c(x, rev(x)), c(y, rep(0, n)),
default.unit="native",
gp=gpar(col=NULL, fill=col))
}
hill()
rdir <- function(n) {
sample(seq(-45, 45, length=10), n)/180*pi
}
grid.text("Long long ago...(ploted by Mr Xia Xiaochao)",
x=.2, y=.51, just="bottom",
default.unit="native",
gp=gpar(col=7, fontface="italic", fontsize=10))
popViewport()
grid.rect()
#######################################################################
par(mfrow=c(2, 2),col=2,bg="light green")
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)
# Don't draw the grid lines : border = NA
par(mar=rep(0, 4))
persp(x, y, z, theta = 135, phi = 30, col = "light blue", scale = FALSE,
ltheta = -120, shade = 0.75, border = NA, box = FALSE)
mtext("persp()", side=3, line=-2)
par(mar=c(3, 3, 2, 0.5))
# Note that example(trees) shows more sensible plots!
N <- nrow(trees)
attach(trees)
# Girth is diameter in inches
symbols(Height, Volume, circles=Girth/24, inches=FALSE,
main="", xlab="", ylab="", bg=grey(Girth/max(Girth)))
mtext("symbols()", side=3, line=0.5)
par(mar=rep(0.5, 4),col.main=4)
contour(x, y, z, asp=1, labcex=0.35, axes=FALSE)
rect(0, 0, 870, 620)
mtext("contour()", side=3, line=-1.5)
image(x, y, z, asp=1, col=grey(0.5 + 1:12/24), xlab="", ylab="", axes=FALSE)
rect(min(x)-5, min(y)-5, max(x)+5, max(y)+5)
mtext("image()", side=3, line=-1.5)
##################################################################################################
#
# Comment:
#
# A bit of mucking around is required to get the second (whole-world)
# map positioned correctly; this provides an example of calling a
# plotting function to perform calculations but do no drawing (see the
# second call to the map() function).
#
# Makes use of the "maps" and "mapproj" packages to draw the maps.
#
library(maps)
par(mar=rep(0, 4))
map("nz", fill=TRUE, col="grey80")
points(174.75, -36.87, pch=16, cex=2)
arrows(172, -36.87, 174, -36.87, lwd=3)
text(172, -36.87, "Auckland ", adj=1, cex=2)
# mini world map as guide
maplocs <- map(projection="sp_mercator", wrap=TRUE, lwd=0.1,
col="grey", ylim=c(-60, 75),
interior=FALSE, orientation=c(90, 180, 0), add=TRUE,
plot=FALSE)
xrange <- range(maplocs$x, na.rm=TRUE)
yrange <- range(maplocs$y, na.rm=TRUE)
aspect <- abs(diff(yrange))/abs(diff(xrange))
# customised to 6.5 by 4.5 figure size
par(fig=c(0.99 - 0.5, 0.99, 0.01, 0.01 + 0.5*aspect*4.5/6.5),
mar=rep(0, 4), new=TRUE)
plot.new()
plot.window(xlim=xrange,
ylim=yrange)
map(projection="sp_mercator", wrap=TRUE, lwd=0.1, ylim=c(-60, 75),
interior=FALSE, orientation=c(90, 180, 0), add=TRUE)
symbols(-.13, -0.8, circles=1, inches=0.1, add=TRUE)
##################################################################################################
#
# Comment:
#
# Code by Arden Miller (Department of Statistics, The University of Auckland).
#
# Lots of coordinate transformations being done "by hand".
# This code is not really reusable; just a demonstration that very
# pretty results are possible if you're sufficiently keen.
#
par(mfrow=c(2, 1), pty="s", mar=rep(1, 4))
# Create plotting region and plot outer circle
plot(c(-1.1, 1.2), c(-1.1, 1.2),
type="n", xlab="", ylab="",
xaxt="n", yaxt="n", cex.lab=2.5)
angs <- seq(0, 2*pi, length=500)
XX <- sin(angs)
YY <- cos(angs)
lines(XX, YY, type="l")
# Set constants
phi1 <- pi*2/9
k1 <- sin(phi1)
k2 <- cos(phi1)
# Create grey regions
obsphi <- pi/12
lambdas <- seq(-pi, pi, length=500)
xx <- cos(pi/2 - obsphi)*sin(lambdas)
yy <- k2*sin(pi/2 - obsphi)-k1 * cos(pi/2 - obsphi)*cos(lambdas)
polygon(xx, yy, col="grey")
lines(xx, yy, lwd=2)
theta1sA <- seq(-obsphi, obsphi, length=500)
theta2sA <- acos(cos(obsphi)/cos(theta1sA))
theta1sB <- seq(obsphi, -obsphi, length=500)
theta2sB <- -acos(cos(obsphi)/cos(theta1sB))
theta1s <- c(theta1sA, theta1sB)
theta2s <- c(theta2sA, theta2sB)
xx <- cos(theta1s)*sin(theta2s+pi/4)
yy <- k2*sin(theta1s)-k1*cos(theta1s)*cos(theta2s+pi/4)
polygon(xx, yy, col="grey")
lines(xx, yy, lwd=2)
xx <- cos(theta1s)*sin(theta2s-pi/4)
yy <- k2*sin(theta1s)-k1*cos(theta1s)*cos(theta2s-pi/4)
polygon(xx, yy, col="grey")
lines(xx, yy, lwd=2)
# Plot longitudes
vals <- seq(0, 7, 1)*pi/8
for(lambda in vals){
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, pi+phi, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=.5)
}
# Grey out polar cap
phi <- 5.6*pi/12
lambdas <- seq(-pi, pi, length=500)
xx <- cos(phi)*sin(lambdas)
yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
polygon(xx, yy, col="grey")
# Plot Latitudes
vals2 <- seq(-2.8, 5.6, 1.4)*pi/12
for(phi in vals2){
if (k1*sin(phi) > k2 * cos(phi))
crit <- pi
else
crit <- acos((-k1*sin(phi))/(k2*cos(phi)))
lambdas <- seq(-crit, crit, length=500)
xx <- cos(phi)*sin(lambdas)
yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
lines(xx, yy, lwd=.5)
}
# Plots axes and label
lines(c(0.00, 0.00), c(k2*sin(pi/2), 1.11), lwd=4)
lines(c(0.00, 0.00), c(-1, -1.12), lwd=4)
a2x <- sin(-pi/4)
a2y <- cos(-pi/4)*(-k1)
lines(c(a2x, 1.5*a2x), c(a2y, 1.5*a2y), lwd=4)
k <- sqrt(a2x^2+a2y^2)
lines(c(-a2x/k, 1.2*(-a2x/k)), c(-a2y/k, 1.2*(-a2y/k)), lwd=4)
a3x <- sin(pi/4)
a3y <- cos(pi/4)*(-k1)
lines(c(a3x, 1.5*a3x), c(a3y, 1.5*a3y), lwd=4)
k <- sqrt(a3x^2+a3y^2)
lines(c(-a3x/k, 1.2*(-a3x/k)), c(-a3y/k, 1.2*(-a3y/k)), lwd=4)
text(0.1, 1.12, expression(bold(X[1])))
text(-1.07, -.85, expression(bold(X[2])))
text(1.11, -.85, expression(bold(X[3])))
# set plot region and draw outer circle
plot(c(-1.1, 1.2), c(-1.1, 1.2),
type="n", xlab="", ylab="",
xaxt="n", yaxt="n", cex.lab=2.5)
angs <- seq(0, 2*pi, length=500)
XX <- sin(angs)
YY <- cos(angs)
lines(XX, YY, type="l")
# set constants
phi1 <- pi*2/9
k1 <- sin(phi1)
k2 <- cos(phi1)
obsphi <- pi/24
# create X2X3 grey region and plot boundary
crit <- acos((-k1*sin(obsphi))/(k2 * cos(obsphi)))
lambdas <- seq(-crit, crit, length=500)
xx1 <- cos(obsphi)*sin(lambdas)
yy1 <- k2*sin(obsphi)-k1 * cos(obsphi)*cos(lambdas)
obsphi <- -pi/24
crit <- acos((-k1*sin(obsphi))/(k2 * cos(obsphi)))
lambdas <- seq(crit, -crit, length=500)
xx3 <- cos(obsphi)*sin(lambdas)
yy3 <- k2*sin(obsphi)-k1 * cos(obsphi)*cos(lambdas)
ang1 <- atan(xx1[500]/yy1[500])
ang2 <- pi+atan(xx3[1]/yy3[1])
angs <- seq(ang1, ang2, length=50)
xx2 <- sin(angs)
yy2 <- cos(angs)
ang4 <- atan(xx1[1]/yy1[1])
ang3 <- -pi+ atan(xx3[500]/yy3[500])
angs <- seq(ang3, ang4, length=50)
xx4 <- sin(angs)
yy4 <- cos(angs)
xxA <- c(xx1, xx2, xx3, xx4)
yyA <- c(yy1, yy2, yy3, yy4)
polygon(xxA, yyA, border="grey", col="grey")
xx1A <- xx1
yy1A <- yy1
xx3A <- xx3
yy3A <- yy3
# create X1X3 grey region and plot boundary
obsphi <- pi/24
crit <- pi/2-obsphi
theta1sA <- c(seq(-crit, crit/2, length=200), seq(crit/2, crit, length=500))
theta2sA <- asin(cos(crit)/cos(theta1sA))
theta1sB <- seq(crit, crit/2, length=500)
theta2sB <- pi-asin(cos(crit)/cos(theta1sB))
theta1s <- c(theta1sA, theta1sB)
theta2s <- c(theta2sA, theta2sB)
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s+pi/4)
xx1 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]+pi/4)
yy1 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]+pi/4)
theta2s <- -theta2s
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s+pi/4)
xx3 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]+pi/4)
yy3 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]+pi/4)
rev <- seq(length(xx3), 1, -1)
xx3 <- xx3[rev]
yy3 <- yy3[rev]
ang1 <- pi+atan(xx1[length(xx1)]/yy1[length(yy1)])
ang2 <- pi+atan(xx3[1]/yy3[1])
angs <- seq(ang1, ang2, length=50)
xx2 <- sin(angs)
yy2 <- cos(angs)
ang4 <- pi+atan(xx1[1]/yy1[1])
ang3 <- pi+atan(xx3[length(xx3)]/yy3[length(yy3)])
angs <- seq(ang3, ang4, length=50)
xx4 <- sin(angs)
yy4 <- cos(angs)
xxB <- c(xx1, -xx2, xx3, xx4)
yyB <- c(yy1, -yy2, yy3, yy4)
polygon(xxB, yyB, border="grey", col="grey")
xx1B <- xx1
yy1B <- yy1
xx3B <- xx3
yy3B <- yy3
# create X1X2 grey region and plot boundary
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s-pi/4)
xx1 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]-pi/4)
yy1 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]-pi/4)
theta2s <- -theta2s
vals <- k1*sin(theta1s)+k2*cos(theta1s)*cos(theta2s-pi/4)
xx3 <- cos(theta1s[vals>=0])*sin(theta2s[vals>=0]-pi/4)
yy3 <- k2*sin(theta1s[vals>=0])-k1*cos(theta1s[vals>=0])*cos(theta2s[vals>=0]-pi/4)
rev <- seq(length(xx3), 1, -1)
xx3 <- xx3[rev]
yy3 <- yy3[rev]
ang1 <- pi+atan(xx1[length(xx1)]/yy1[length(yy1)])
ang2 <- pi+atan(xx3[1]/yy3[1])
angs <- seq(ang1, ang2, length=50)
xx2 <- sin(angs)
yy2 <- cos(angs)
ang4 <- pi+atan(xx1[1]/yy1[1])
ang3 <- pi+atan(xx3[length(xx3)]/yy3[length(yy3)])
angs <- seq(ang3, ang4, length=50)
xx4 <- sin(angs)
yy4 <- cos(angs)
xx <- c(xx1, -xx2, xx3, xx4)
yy <- c(yy1, -yy2, yy3, yy4)
polygon(xx, yy, border="grey", col="grey")
xx1C <- xx1
yy1C <- yy1
xx3C <- xx3
yy3C <- yy3
# plot boundaries to grey regions
lines(xx1C[2:45], yy1C[2:45], lwd=2)
lines(xx1C[69:583], yy1C[69:583], lwd=2)
lines(xx1C[660:1080], yy1C[660:1080], lwd=2)
lines(xx3C[13:455], yy3C[13:455], lwd=2)
lines(xx3C[538:1055], yy3C[538:1055], lwd=2)
lines(xx3C[1079:1135], yy3C[1079:1135], lwd=2)
lines(xx1A[6:113], yy1A[6:113], lwd=2)
lines(xx1A[153:346], yy1A[153:346], lwd=2)
lines(xx1A[389:484], yy1A[389:484], lwd=2)
lines(xx3A[1:93], yy3A[1:93], lwd=2)
lines(xx3A[140:362], yy3A[140:362], lwd=2)
lines(xx3A[408:497], yy3A[408:497], lwd=2)
lines(xx1B[2:45], yy1B[2:45], lwd=2)
lines(xx1B[69:583], yy1B[69:583], lwd=2)
lines(xx1B[660:1080], yy1B[660:1080], lwd=2)
lines(xx3B[13:455], yy3B[13:455], lwd=2)
lines(xx3B[538:1055], yy3B[538:1055], lwd=2)
lines(xx3B[1079:1135], yy3B[1079:1135], lwd=2)
# Plot longitudes
vals <- seq(-7, 8, 1)*pi/8
for(lambda in vals){
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, 5.6*pi/12, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=.5)
}
# Plot Latitudes
# vals2 <- seq(-2.8, 5.6, 1.4)*pi/12
vals2 <- c(-1.5, 0, 1.5, 3.0, 4.5, 5.6)*pi/12
for(phi in vals2){
if (k1*sin(phi) > k2 * cos(phi))
crit <- pi
else
crit <- acos((-k1*sin(phi))/(k2*cos(phi)))
lambdas <- seq(-crit, crit, length=500)
xx <- cos(phi)*sin(lambdas)
yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
lines(xx, yy, lwd=.5)
}
# create lines for X1X2- and X1X3-planes
lambda <- pi/4
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, pi+phi, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=2)
lambda <- 3*pi/4
sl <- sin(lambda)
cl <- cos(lambda)
phi <- atan(((0-1)*k2*cl)/(k1))
angs <- seq(phi, pi+phi, length=500)
xx <- cos(angs)*sl
yy <- k2*sin(angs)-k1*cos(angs)*cl
lines(xx, yy, lwd=2)
# create line for X2X3-plane
phi <- 0
crit <- acos((-k1*sin(phi))/(k2 * cos(phi)))
lambdas <- seq(-crit, crit, length=500)
xx <- cos(phi)*sin(lambdas)
yy <- k2*sin(phi)-k1 * cos(phi)*cos(lambdas)
lines(xx, yy, lwd=2)
# create axes
lines(c(0.00, 0.00), c(k2*sin(pi/2), 1.11), lwd=4)
lines(c(0.00, 0.00), c(-1, -1.12), lwd=4)
a2x <- sin(-pi/4)
a2y <- cos(-pi/4)*(-k1)
lines(c(a2x, 1.5*a2x), c(a2y, 1.5*a2y), lwd=4)
a3x <- sin(pi/4)
a3y <- cos(pi/4)*(-k1)
lines(c(a3x, 1.5*a3x), c(a3y, 1.5*a3y), lwd=4)
k <- sqrt(a3x^2+a3y^2)
lines(c(-a3x/k, 1.2*(-a3x/k)), c(-a3y/k, 1.2*(-a3y/k)), lwd=4)
k <- sqrt(a2x^2+a2y^2)
lines(c(-a2x/k, 1.2*(-a2x/k)), c(-a2y/k, 1.2*(-a2y/k)), lwd=4)
# add text
text(-1.07, -.85, expression(bold(X[2])))
text(1.11, -.85, expression(bold(X[3])))
text(0.1, 1.12, expression(bold(X[1])))
lines(XX, YY, type="l")
#########################################################################################################
install.packages("pixmap")
library(pixmap)
moonPhase <- function(x, y, phase, size=.07) {
# phase 1: first quarter
# 2: full
# 3: last quarter
# 4: new
# size is in inches
n <- 17
angle <- seq(0, 2*pi, length=n)
xx <- x + cos(angle)*xinch(size)
yy <- y + sin(angle)*yinch(size)
if (phase == 4)
fill <- "black"
else
fill <- "white"
polygon(xx, yy, col=fill)
if (phase == 1)
polygon(xx[(n/4):(n*3/4) + 1],
yy[(n/4):(n*3/4) + 1],
col="black")
if (phase == 3)
polygon(xx[c(1:(n/4 + 1), (n*3/4 + 1):n)],
yy[c(1:(n/4 + 1), (n*3/4 + 1):n)],
col="black")
}
# Data from Land Information New Zealand
# http://hydro.linz.govt.nz
# + 1 for daylight saving
hours <- c(18, 18, 19, 20, 21, 22, 23,
0, 1, 2, 3, 4, 5, 5,
6, 7, 8, 9, 10, 11, 12,
13, 13, 14, 15, 15, 16) + 1
hours[7] <- 0 # 23 + 1 = 0 on a 24-hour clock
mins <- c(9, 57, 49, 46, 48, 54, 59,
59, 52, 41, 28, 14, 1, 52,
36, 43, 41, 39, 38, 35, 26,
10, 49, 26, 2, 39, 16)
lowTideDate <- ISOdatetime(2005, 2, c(1:6,8:28),
hours, mins, 0)
lowTideHour <- ISOdatetime(2005, 2, 1,
hours, mins, 0)
phases <- ISOdatetime(2005, 2, c(2, 9, 16, 24),
c(19, 10, 12, 16) + 1,
c(28, 30, 16, 55), 0)
mainHours <- ISOdatetime(2005, 2, 1,
c(0, 4, 8, 12, 16, 20, 23),
c(rep(0, 6), 59),
c(rep(0, 6), 59))
# Original image from NASA
# http://grin.hq.nasa.gov/ABSTRACTS/GPN-2000-000473.html
moon <- read.pnm(system.file(file.path("Images",
"GPN-2000-000473halfsize.pnm"),
package="RGraphics"))
par(pty="s", xaxs="i", yaxs="i", cex=.7)
plot.new()
addlogo(moon, 0:1, 0:1, asp=1)
par(new=TRUE, xaxs="r", yaxs="r", las=1)
plot(lowTideDate, lowTideHour, type="n",
ylim=range(mainHours), axes=FALSE, ann=FALSE)
# dashed reference lines
midday <- ISOdatetime(2005, 2, 1, 12, 0, 0)
abline(h=midday, v=phases,
col="white", lty="dashed")
# grey "repeat" of tide info to show gradient
lines(lowTideDate[6:7],
c(ISOdatetime(2005, 1, 31,
hours[6], mins[6], 0),
lowTideHour[7]),
lwd=2, col="grey50")
points(lowTideDate[6:7],
c(ISOdatetime(2005, 1, 31,
hours[6], mins[6], 0),
lowTideHour[7]),
pch=16, cex=.7, col=7) ##grey50
for (subset in list(1:6, 7:27)) {
lines(lowTideDate[subset], lowTideHour[subset],
lwd=2, col="light green")
points(lowTideDate[subset], lowTideHour[subset],
pch=16, cex=.7, col="light blue")
}
box()
axis.POSIXct(1, lowTideDate)
axis.POSIXct(2, at=mainHours, format="%H:%M")
mtext("Time of Low Tide (NZDT)", side=2, line=4, las=0)
mtext("Auckland, New Zealand 2005", side=1, line=3)
axis(3, at=phases, labels=FALSE)
par(xpd=NA,bg=5)
ymax <- par("usr")[4]
for (i in 1:4)
moonPhase(phases[i], ymax + yinch(.2), c(3, 4, 1, 2)[i])
mtext("Phases of the Moon", side=3, line=3)
##################################################################################################
par(mar=rep(0, 4), lwd=0.1,bg=7)
z <- 2 * volcano
x <- 10 * (1:nrow(z))
y <- 10 * (1:ncol(z))
trans <- persp(x, y, z, theta = 135, phi = 30,
scale = FALSE, ltheta = -120,
# shade=0.5, border=NA,
box = FALSE)
box(col="grey", lwd=1)
trans3d <- function(x,y,z,pmat) {
tmat <- cbind(x,y,z,1)%*% pmat
tmat[,1:2] / tmat[,4]
}
summit <- trans3d(x[20], y[31], max(z), trans)
points(summit[1], summit[2], pch=16)
summitlabel <- trans3d(x[20], y[31], max(z) + 50, trans)
text(summitlabel[1], summitlabel[2], "Summit")
drawRoad <- function(x, y, z, trans) {
road <- trans3d(x, y, z, trans)
lines(road[,1], road[,2], lwd=5)
lines(road[,1], road[,2], lwd=3, col="grey")
}
with(volcano.summitRoad,
drawRoad(srx, sry, srz, trans))
with(volcano.upDownRoad,
{
clipudx <- udx
clipudx[udx < 230 & udy < 300 |
udx < 150 & udy > 300] <- NA
drawRoad(clipudx, udy, udz, trans)
})
with(volcano.accessRoad,
drawRoad(arx, ary, arz, trans))
#############################################################################################
par(mfrow=c(3,1),lwd=0.1,bg=5)
grid.rect(gp=gpar(col=2))
suffix <- c("even", "odd")
for (i in 1:8)
grid.circle(name=paste("circle.", suffix[i %% 2 + 1],
sep=" "),
r=(9 - i)/20,
gp=gpar(col=7, fill=grey(i/10)))
grid.edit("circle.odd", gp=gpar(fill="grey10"),
global=TRUE)
grid.edit("circle", gp=gpar(col="grey", fill="grey90"),
grep=TRUE, global=TRUE)
############################################################################################
grid.rect(gp=gpar(col="grey"))
pushViewport(viewport(w=0.8, h=0.8))
offset <- unit(3, "mm")
grid.text("CIE-LUV\nchroma = 40\nluminance = 80",
gp=gpar(fontfamily="mono"))
t <- seq(0, 2*pi, length=7)[-7]
x <- 0.5 + .4*cos(t)
y <- 0.5 + .4*sin(t)
rad <- .13
cols <- hcl(t/pi*180, 40, 80)
grid.circle(x, y, r=unit(rad, "npc"),
gp=gpar(fill=cols))
labels <- c("r", "g", "b")
rgbcols <- col2rgb(cols)
hbars <- seq(50, 250, 50)
for (i in 1:6) {
pushViewport(viewport(x[i], y[i], width=.2, height=.2,
layout=grid.layout(2, 3,
heights=unit(c(1, 1), c("null", "lines")))))
for (j in 1:3) {
pushViewport(viewport(layout.pos.col=j,
layout.pos.row=1))
grid.text(labels[j], y=unit(-.5, "lines"),
gp=gpar(fontfamily="mono"))
grid.lines(x=c(.1, .9), y=0)
grid.rect(y=0, width=.6, height=rgbcols[j, i]/300,
just="bottom",
gp=gpar(fill=NA))
popViewport()
}
popViewport()
}
popViewport()
扫一扫
294
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
