[R] Rose diagrams in R?
David Finlayson
david_p_finlayson at hotmail.com
Mon Nov 26 04:22:52 CET 2001
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list