Dan,
May I ask you to report this as an issue on the oce webpage,
so that others can see the discussion? (The "R help" is
perhaps not the right place to report bugs ... and, yes, this
is a bug.)
http://code.google.com/p/r-oce/issues/list
Yes, I will.
A possible solution is to edit the "windrose.R" file and
replace the first function with what I am putting below. But
I am actually not sure on the "if (theta..." line, and need
to think about that some more. Also, please note that
official releases of OCE take some time (typically a week) so
that's why I always suggest a workaround, first. As I'm sure
you're aware, you could just "source" a file containing the
code below, after doing the "library(oce)", and that will
work until a new release comes out. I have been holding off
on a release for quite a while, since I am working on code to
handle acoustic-doppler data from about a half-dozen
instruments, and I try to avoid letting anyone get caught out
by using code that is still in development.
as.windrose <- function(x, y, dtheta = 15) {
dt <- dtheta * pi / 180
dt2 <- dt / 2
R <- sqrt(x^2 + y^2)
angle <- atan2(y, x)
L <- max(R, na.rm=TRUE)
nt <- 2 * pi / dt
theta <- count <- mean <- vector("numeric", nt)
fivenum <- matrix(0, nt, 5)
for (i in 1:nt) {
theta[i] <- i * dt
if (theta[i] <= pi)
inside <- (angle < (theta[i] + dt2)) &
((theta[i] - dt2) <= angle)
else {
inside <- ((2*pi+angle) < (theta[i] + dt2)) & ((theta[i]
- dt2) <= (2*pi+angle))
}
count[i] <- sum(inside)
mean[i] <- mean(R[inside], na.rm=TRUE)
fivenum[i,] <- fivenum(R[inside], na.rm=TRUE)
}
data <- list(n=length(x), x.mean=mean(x, na.rm=TRUE),
y.mean=mean (y, na.rm=TRUE), theta=theta*180/pi, count=count,
mean=mean,
fivenum=fivenum)
metadata <- list(dtheta=dtheta)
log <- processing.log.item(paste(deparse(match.call()), sep="",
collapse=""))
res <- list(data=data, metadata=metadata, processing.log=log)
class(res) <- c("windrose", "oce")
res
}
I decided to not use the windrose object created by as.windrose().
Instead, I wrote my own function for plotting wind speed and wind
direction data, using some of the features from plot.windrose(). The
data I work with is usually provided in mph or m/s and azimuth
degrees,
so I found it more convenient to provide those directly as arguments
to
the function. Three utility functions come first, followed by the
main
one of interest, PLOT.WINDROSE2().
degrees <- function(radians) {
return(radians * 180 / pi)
}
radians <- function(degrees) {
return(degrees * pi / 180)
}
# Define a function that converts compass headings (bearings,
azimuths)
into geometric angles for trigonometry.
# Examples: For a heading = 0 degrees, vector points due north, and
angle = 90 deg.
# For a heading = 90 degrees, vector points east, and angle = 0 deg.
convert.heading.angle <- function(heading) {
num.heading <- length(heading)
angles <- rep(NA, num.heading)
for(i in 1:num.heading) {
angles[i] <-
ifelse((heading[i] >=0 && heading[i] <= 90), 90 - heading[i],
ifelse((heading[i] >=90 && heading[i] <= 180), 360 -
(heading[i] - 90),
ifelse((heading[i] >=180 && heading[i] <= 270), 180
+ (270 - heading[i]), 90 + (360 - heading[i]))
)
)
}
return(angles)
} # end of convert.heading.angle()
# A function to plot windroses, with similar output capabilities as
that
of plot.windrose() function
# in the oce package. This function has not been extensively tested.
# Theta should be in ascending order, starting at some azimuth > 0
degrees.
# The oce package was written by Dan Kelley.
# Arguments:
# 1) vector of windspeeds
# 2) vector of wind directions (azimuths, degrees; 0=calm, 45=NE,
90=E,
180=S, 270=W, 360=N)
# AZIMUTH DEGREES ARE COMPASS BEARINGS
# 3) vector of windrose sector centers (azimuths, degrees)
# 4) type of windrose plot: count, mean, median, or fivenum
# 5) color vector for use as in plot.windrose()
# 6) plot title
PLOT.WINDROSE2 <- function(ws, az, theta, type=c("count", "mean",
"median", "fivenum"), col, title, ...) {
type <- match.arg(type)
nt <- length(theta) # number of sectors
az.ar <- radians(az) # convert azimuth degrees to radians
t <- radians(theta) # convert sector centers to radians
dt <- t[2] - t[1] # angle range for each sector (radians)
if(dt < 0) stop("Center angles for theta must be in increasing
order.\n")
if((sum(diff(theta)) + theta[2]-theta[1]) != 360) stop("Sectors do
not
add up to 360 degrees.\n")
dt2 <- dt / 2 # radians
sector.boundaries <- c(t - dt2, t[nt] + dt2) # sector boundaries in
azimuth radians
# bin the azimuths into the sectors
ind.sector <- ifelse(az.ar < sector.boundaries[1], nt,
findInterval(az.ar, sector.boundaries))
# Bin the azimuth data into the defined sectors
calm <- ifelse(az == 0, T, F) # logical vector to indicate whether
calm or not
# For each sector, compute the count and mean, and the components of
fivenum:
# min, 1st quartile, median, 3rd quartile, and max.
x <- list() # make a list like plot.windrose() uses
x$data <- list(count=integer(nt), mean=numeric(nt),
median=numeric(nt), fivenum=matrix(nrow=nt, ncol=5))
for(i in 1:nt) {
ind <- which(ind.sector == i & !calm)
x$data$count[i] <- length(ind)
x$data$mean[i] <- mean(ws[ind])
x$data$fivenum[i,] <- fivenum(ws[ind])
}
# Plot setup
plot.new()
xlim <- ylim <- c(-1, 1)
par(xpd=NA)
plot.window(xlim, ylim, "", asp = 1)
if (missing(col)) {
col <- c("red", "pink", "blue", "darkgray")
} else {
if (length(col) != 4) stop("'col' should be a list of 4 colors")
}
# Draw circle and radii
tt <- seq(0, 2*pi, length.out=100)
px <- cos(tt)
py <- sin(tt)
lines(px, py, col=col[4]) # draw the circle
# Convert sector boundaries from azimuth radians to geometric radians
for drawing
t.deg <- degrees(t)
t <- radians(convert.heading.angle(t.deg))
for (i in 1:nt) {
lines(c(0, cos(t[i] - dt2)), c(0, sin(t[i] - dt2)), lwd=0.5,
col=col[4]) # draw the radius
}
text( 0, -1, "S", pos=1)
text(-1, 0, "W", pos=2)
text( 0, 1, "N", pos=3)
text( 1, 0, "E", pos=4)
# Draw rose in a given type
if (type == "count") {
max <- max(x$data$count, na.rm=TRUE)
for (i in 1:nt) {
r <- x$data$count[i] / max
xlist <- c(0, r * cos(t[i] - dt2), r * cos(t[i] + dt2), 0)
ylist <- c(0, r * sin(t[i] - dt2), r * sin(t[i] + dt2), 0)
polygon(xlist, ylist,col=col[1], border=col[1])
}
title(title)
} else if (type == "mean") {
max <- max(x$data$mean, na.rm=TRUE)
for (i in 1:nt) {
r <- x$data$mean[i] / max
xlist <- c(0, r * cos(t[i] - dt2), r * cos(t[i] + dt2), 0)
ylist <- c(0, r * sin(t[i] - dt2), r * sin(t[i] + dt2), 0)
polygon(xlist, ylist,col=col[1], border=col[1])
}
title(title)
} else if (type == "median") {
max <- max(x$data$fivenum[,3], na.rm=TRUE)
for (i in 1:nt) {
r <- x$data$fivenum[i,3] / max
xlist <- c(0, r * cos(t[i] - dt2), r * cos(t[i] + dt2), 0)
ylist <- c(0, r * sin(t[i] - dt2), r * sin(t[i] + dt2), 0)
polygon(xlist, ylist,col=col[1], border=col[1])
}
title(title)
} else if (type == "fivenum") {
max <- max(x$data$fivenum[,5], na.rm=TRUE)
for (i in 1:nt) {
for (j in 2:5) {
tm <- t[i] - dt2
tp <- t[i] + dt2
r0 <- x$data$fivenum[i, j-1] / max
r <- x$data$fivenum[i, j ] / max
xlist <- c(r0 * cos(tm), r * cos(tm), r * cos(tp), r0 *
cos(tp))
ylist <- c(r0 * sin(tm), r * sin(tm), r * sin(tp), r0 *
sin(tp))
thiscol <- col[c(2,1,1,2)][j-1]
polygon(xlist, ylist, col=thiscol, border=col[4])
}
r <- x$data$fivenum[i, 3] / max
lines(c(r * cos(tm), r * cos(tp)), c(r * sin(tm), r *
sin(tp)), col="blue", lwd=2)
}
title(title)
}
invisible()
} # end of PLOT.WINDROSE2()
# Example:
x11()
az <- c(1, 22.5, 45, 90, 180, 270, 315, 359, 360, 192, 135)
ws <- rep(1, length(az))
theta <- seq(22.5, 360, by=22.5)
title <- "Test"
PLOT.WINDROSE2(ws, az, theta, type="mean", title=title)
--Scott Waichler
scott.waich...@pnl.gov
