[R] Output from as.windrose() in oce package baffles me
Waichler, Scott R
Scott.Waichler at pnl.gov
Sat Sep 5 00:06:40 CEST 2009
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.waichler at pnl.gov
More information about the R-help
mailing list