[R] Rose diagrams in R?

[David Finlayson]

Thanks all for the help.  I was inspired to write my own version.  It's a
hybrid of Joerg's, Ben's and my own doodling around.  It could still use
some finishing touches (like rings don't make sense when you have no numeric
scale against which to compare them) but it will function for my purposes.
It differs from Ben and Joerg's mainly in that it accepts the normal
graphics calls, only labels the top 10 percent of "petals" and it adds
cardinal direction labels.  Also, it rotates the plot so that 000 is
pointing towards the top of the page (normal convention in the West for map
data).  I don't believe that Ben's did this as written.  Finally, it makes
an effort to scale things correctly and returns a histogram object.

I was impressed by how (relatively) easy it was to implement.

David Finlayson

Code Follows:

# Directional Vector Histogram (Rose Diagrams).
# -------------------------------------------
# David Finlayson (with help from Joerg Maeder and Ben Bolker)
# david_p_finlayson at hotmail.com
# Version 1.0
# November 23, 2001
#
# Use for plotting directional data such as wind direction or the
# angles of imbricated pebbles in rivers and streams.  This is basically
# an extension of the hist(x) function though I did not implement all of
# hist.  I have placed limits on the range of bins so that they always
# fall within 0 and 360 (i.e. directions of the compass).  The standard
# color and line adjustment commands work as well but you will need to add
# annotation (i.e.. main, xlab, ylab) separately (see par).
#
# bins: Approximate number of bins see hist() function for details
# rscale: Ring Scale, the approximate number of rings for scaling see
pretty()
#            NULL value will call pretty() with default number of rings
# labels:  (T/F) draw labels for the top 10% largest petals and the cardinal
directions?
# rings: (T/F) draw scale rings?
#
# example:
# test <- runif(30) * 360
# rose(test)
# rose(test, bins=10, rscale=2, labels=TRUE, rings=TRUE, col="cyan", lwd=2)
#

rose <- function(x, bins=36, rscale=NULL, labels=TRUE, rings=TRUE, ...){

 ### Ensure that this is directional data (0-360)
 if (max(x) > 360 || min(x) < 0) {
  stop("Data is out of range (0 - 360)")
 }

 ### Histogram Data
 histogram.out <- hist(x, breaks=seq(0,360,length=bins+1), plot=FALSE)
 pieMid <- histogram.out$mids  # mid points of bins
 pieCount <- histogram.out$counts # count in each bin
 pieWidth <- 360 / length(pieMid) # width of each bin
 pieFreq <- histogram.out$density * pieWidth

 ### Initialize Plot Are
 oldpar <- par()
 par(pty="s")

 if (!is.null(rscale)){
  rscale <- pretty(pieFreq, rscale)
   } else {
  rscale <- pretty(pieFreq)
 }

 plotLimits <- c(-max(rscale), max(rscale)) * 1.04

 plot(0,0,
  ylim=plotLimits,
  xlim=plotLimits,
  axes=FALSE,
  xlab="",
  ylab="")

 abline(h=0)
 abline(v=0)

 ### Plot Rings
 if (rings == TRUE) {
  symbols(rep(0,length(rscale)), rep(0,length(rscale)),
   circles=rscale,inches=FALSE,add=TRUE)
 }

 ### Plot Rose

 # Helper function to draw the pie slices
 pie <- function(h, k, direction, spread, magnitude, ...){

  # adjust coordinate from compass to polar
  direction <- 360 - direction + 90

  # calculate theta start and stop
  start <- (direction + 0.5 * spread) * pi / 180
  stop  <- (direction - 0.5 * spread) * pi / 180

  # vertices
  x1 <- h
  y1 <- k

  x2 <- (magnitude * cos(start)) + h
  y2 <- (magnitude * sin(start)) + k

  x3 <- (magnitude * cos(stop)) + h
  y3 <- (magnitude * sin(stop)) + k

  # build pie slice
  x <- c(x1, x2, x3, x1)
  y <- c(y1, y2, y3, y1)

  polygon(x,y, ...)
 }

 # plot the slices
 for (i in 1:length(pieMid)){
   pie(0, 0, pieMid[i], pieWidth, pieFreq[i], ...)
 }

 # Plot Axes Labels

 if (labels == TRUE) {

     mtext("N", side=3, line=1)
     mtext("E", side=4, line=1)
     mtext("S", side=1, line=1)
     mtext("W", side=2, line=1)

     ### Plot top frequency labels

     # calculate indices of top 10 percent of bins
     pie10percent <- round(length(pieMid) * 0.1) + 1 # how many is 10
percent
     pieRank <- length(pieCount) - rank(pieCount)   # rank bins
     top10 <- which( pieRank < pie10percent )  # index to top 10 percent

     # Plot labels for these bins
     theta <- 360 - pieMid[top10] + 90 # compass to polar angles
     theta <- theta * pi/180  # degrees to rads

     x <- pieFreq[top10] * cos(theta)
     y <- pieFreq[top10] * sin(theta)

     text(x, y, format(pieFreq[top10], digits=2))

     ### Reset the par
     par <- oldpar

     ### Return Histogram Object
     histogram.out
 }
}


----- Original Message -----
From: "Paul Murrell" <p.murrell at auckland.ac.nz>
To: "Joerg Maeder" <joerg.maeder at ethz.ch>; "David Finlayson"
<david_p_finlayson at hotmail.com>; "'R-help at lists.R-project.org'"
<R-help at stat.math.ethz.ch>
Sent: Sunday, November 25, 2001 11:28 AM
Subject: Re: [R] Rose diagrams in R?


> Hi
>
>
> > i wrote a small function (rose) for that problem. It accept datas
> > between 0 and 360 and draw something
> > like you want. You also can
> >
> > rose <- function(data,step=30,main='wind rose'){
> >   deg2rad <- 180/pi
> >   data <- (data+step/2)%%360# Values like 359 go to Sector 0
> >   histdata <- hist(data,breaks=seq(0,360,by=step),plot=F) #use hist for
> > counting
> >   counts <- histdata$counts
> >   maxcount <- max(counts)
> >   mids <- (histdata$mids-step/2)/deg2rad
> >   step <- step/deg2rad
> >   plot(c(-1,1),c(-1,1),,xlab='',ylab='',
> >        main=main,xaxt='n',yaxt='n',pch=' ')
> >   for (i in 1:length(counts)){
> >     w1 <- mids[i]-step/2
> >     w2 <- mids[i]+step/2
> >     lines(counts[i]/maxcount*c(0,sin(w1),sin(w2),0),
> >           counts[i]/maxcount*c(0,cos(w1),cos(w2),0))#draw sector
> >     text(sin(mids[i]),cos(mids[i]),counts[i])
> >   }
> >   names(counts) <- round(mids*deg2rad,3)
> >   counts
> > }
> >
> > #Test with 500 values between 0 and 360 (uniform distribution)
> > rose(runif(500)*360,360/16)
>
>
> You might want to include a par(pty="s") in there somewhere (or something
> equivalent) to make sure that the plot comes out "square".
>
> Paul
>
>
>

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._