[R] Output from as.windrose() in oce package baffles me

[Dan Kelley]

Thanks for posting your code, Scott.  I'll have a look at the windrose  
in OCE sometime next week, to see if there is anything I can do that  
would have made your task easier.

I've never used windroses in my own work, so that may explain why the  
interface is clunky!  Also, it's not completely clear to me what plots  
"should" look like.  This relates a bit to the issue of pie charts  
being difficult to interpret by eye, for one thing.  One thing that I  
*have* seen in my field is a polar diagram of wind speed, typically  
drawn as an image.  But, there again, there is the issue that radial  
lines diverge, which confuses the eye somewhat.  The area of a  
coloured patch of wind of a patch gets larger as the radius, or wind  
speed, increases.   Now, that may make sense, if you're interested in  
wind stress (proportional to the square of wind).

The bottom line is that it seems to me that the binning into delta- 
theta classes is problematic.  My tendency would be to divide the wind  
into dx and dy classes, and to plot a two-dimensional histogram.  Of  
course, that won't look at all like a wind-rose!

Thanks again.  Dan.

On 2009-09-04, at 7:06 PM, Waichler, Scott R wrote:

> 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
>
>