R help - Sep 2009 - Output from as.windrose() in oce package baffles me

I'm having trouble understanding the output from as.windrose().  For one
thing, data on a boundary between sectors seem to be left out of the
counts.  I assume that explains the missing point in the output below
(angle 45).  Shouldn't one side of each sector interval be open, to
include values such as my 45 in the example?  Also, why does the angle
180 in my input apparently not result in a count of 1 for the output for
180?  The $data$theta value of 180 refers to a sector centered at 180
whose angle is 90, correct?

library(oce)

# azimuths <- c(22.5, 45, 90, 180, 270)  # 0 and 360=north, 90=east,
180=south, 270=west
angles <- c(67.5, 45, 0, 270, 180)  # geometric equivalents of azimuths
radians <- angles * pi /180
x <- round(cos(radians), 6)
y <- round(sin(radians), 6)
w <- as.windrose(x, y, dtheta = 90) # make a windrose object, oce
package

# part of windrose object w:
$data$theta
[1]  90 180 270 360

$data$count
[1] 1 0 1 1

Thanks,
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler at pnl.gov

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