R/windRose.R
In davidcarslaw/ggopenair: Tools for the Analysis of Air Pollution Data
pollutionRose <- function(mydata, pollutant = "nox", key.footer = pollutant,
key.position = "right", key = TRUE,
breaks = 6, paddle = FALSE, seg = 0.9, normalise = FALSE,
...)
{
## extra args setup
extra <- list(...)
## check to see if two met data sets are being compared.
## if so, set pollutant to one of the names
if ("ws2" %in% names(extra)) {
pollutant <- extra$ws
if (missing(breaks)) breaks <- NA
}
if (is.null(breaks)) breaks <- 6
if (is.numeric(breaks) & length(breaks) == 1) {
## breaks from the minimum to 90th percentile, which generally gives sensible
## spacing for skewed data. Maximum is added later.
breaks <- unique(pretty(c(min(mydata[[pollutant]], na.rm = TRUE),
quantile(mydata[[pollutant]], probs = 0.9, na.rm = TRUE),
breaks)))
}
windRose(mydata, pollutant = pollutant, paddle = paddle, seg = seg,
key.position = key.position, key.footer = key.footer, key = key,
breaks = breaks, normalise = normalise, ...)
}
##' Traditional wind rose plot and pollution rose variation
##'
##' The traditional wind rose plot that plots wind speed and wind direction by
##' different intervals. The pollution rose applies the same plot structure but
##' substitutes other measurements, most commonly a pollutant time series, for
##' wind speed.
##'
##' For \code{windRose} data are summarised by direction, typically by 45 or 30
##' (or 10) degrees and by different wind speed categories. Typically, wind
##' speeds are represented by different width "paddles". The plots show the
##' proportion (here represented as a percentage) of time that the wind is from
##' a certain angle and wind speed range.
##'
##' By default \code{windRose} will plot a windRose in using "paddle" style
##' segments and placing the scale key below the plot.
##'
##' The argument \code{pollutant} uses the same plotting structure but
##' substitutes another data series, defined by \code{pollutant}, for wind
##' speed.
##'
##' The option \code{statistic = "prop.mean"} provides a measure of the
##' relative contribution of each bin to the panel mean, and is intended for
##' use with \code{pollutionRose}.
##'
##' \code{pollutionRose} is a \code{windRose} wrapper which brings
##' \code{pollutant} forward in the argument list, and attempts to sensibly
##' rescale break points based on the \code{pollutant} data range by by-passing
##' \code{ws.int}.
##'
##' By default, \code{pollutionRose} will plot a pollution rose of \code{nox}
##' using "wedge" style segments and placing the scale key to the right of the
##' plot.
##'
##' It is possible to compare two wind speed-direction data sets using
##' \code{pollutionRose}. There are many reasons for doing so e.g. to
##' see how one site compares with another or for meteorological model
##' evaluation. In this case, \code{ws} and \code{wd} are considered
##' to the the reference data sets with which a second set of wind
##' speed and wind directions are to be compared (\code{ws2} and
##' \code{wd2}). The first set of values is subtracted from the second
##' and the differences compared. If for example, \code{wd2} was
##' biased positive compared with \code{wd} then \code{pollutionRose}
##' will show the bias in polar coordinates. In its default use, wind
##' direction bias is colour-coded to show negative bias in one colour
##' and positive bias in another.
##'
##' @usage windRose(mydata, ws = "ws", wd = "wd", ws2 = NA, wd2 = NA,
##' ws.int = 2, angle = 30, type = "default", bias.corr = TRUE, cols = "default",
##' grid.line = NULL, width = 1, seg = NULL, auto.text = TRUE, breaks
##' = 4, offset = 10, normalise = FALSE, max.freq = NULL, paddle = TRUE, key.header =
##' NULL, key.footer = "(m/s)", key.position = "bottom", key = TRUE,
##' dig.lab = 5, statistic = "prop.count", pollutant = NULL, annotate
##' = TRUE, angle.scale = 315, border = NA, ...)
##'
##'
##' pollutionRose(mydata, pollutant = "nox", key.footer = pollutant,
##' key.position = "right", key = TRUE, breaks = 6, paddle = FALSE,
##' seg = 0.9, normalise = FALSE, ...)
##'
##'
##' @aliases windRose pollutionRose
##' @param mydata A data frame containing fields \code{ws} and \code{wd}
##' @param ws Name of the column representing wind speed.
##' @param wd Name of the column representing wind direction.
##' @param ws2 The user can supply a second set of wind speed and wind
##' direction values with which the first can be compared. See details
##' below for full explanation.
##' @param wd2 see \code{ws2}.
##' @param ws.int The Wind speed interval. Default is 2 m/s but for low met
##' masts with low mean wind speeds a value of 1 or 0.5 m/s may be better.
##' Note, this argument is superseded in \code{pollutionRose}. See
##' \code{breaks} below.
##' @param angle Default angle of \dQuote{spokes} is 30. Other potentially useful
##' angles are 45 and 10. Note that the width of the wind speed interval may
##' need adjusting using \code{width}.
##' @param type \code{type} determines how the data are split
##' i.e. conditioned, and then plotted. The default is will produce a
##' single plot using the entire data. Type can be one of the built-in
##' types as detailed in \code{cutData} e.g. \dQuote{season},
##' \dQuote{year}, \dQuote{weekday} and so on. For example, \code{type
##' = "season"} will produce four plots --- one for each season.
##'
##' It is also possible to choose \code{type} as another variable in the data
##' frame. If that variable is numeric, then the data will be split into four
##' quantiles (if possible) and labelled accordingly. If type is an existing
##' character or factor variable, then those categories/levels will be used
##' directly. This offers great flexibility for understanding the variation
##' of different variables and how they depend on one another.
##'
##' Type can be up length two e.g. \code{type = c("season", "weekday")} will
##' produce a 2x2 plot split by season and day of the week. Note, when two
##' types are provided the first forms the columns and the second the rows.
##' @param bias.corr When \code{angle} does not divide exactly into
##' 360 a bias is introduced in the frequencies when the wind
##' direction is already supplied rounded to the nearest 10 degrees,
##' as is often the case. For example, if \code{angle = 22.5}, N, E,
##' S, W will include 3 wind sectors and all other angles will be
##' two. A bias correction can made to correct for this problem. A
##' simple method according to Applequist (2012) is used to adjust the
##' frequencies.
##' @param cols Colours to be used for plotting. Options include
##' \dQuote{default}, \dQuote{increment}, \dQuote{heat}, \dQuote{jet},
##' \dQuote{hue} and user defined. For user defined the user can
##' supply a list of colour names recognised by R (type
##' \code{colours()} to see the full list). An example would be
##' \code{cols = c("yellow", "green", "blue", "black")}.
##' @param grid.line Grid line interval to use. If \code{NULL}, as in default,
##' this is assigned by \code{windRose} based on the available data range.
##' However, it can also be forced to a specific value, e.g.
##' \code{grid.line = 10}.
##' @param width For \code{paddle = TRUE}, the adjustment factor for width of
##' wind speed intervals. For example, \code{width = 1.5} will make the
##' paddle width 1.5 times wider.
##' @param seg For \code{pollutionRose} \code{seg} determines with
##' width of the segments. For example, \code{seg = 0.5} will produce
##' segments 0.5 * \code{angle}.
##' @param auto.text Either \code{TRUE} (default) or \code{FALSE}. If
##' \code{TRUE} titles and axis labels will automatically try and
##' format pollutant names and units properly e.g. by subscripting
##' the \sQuote{2} in NO2.
##' @param breaks Most commonly, the number of break points for wind
##' speed in \code{windRose} or pollutant in \code{pollutionRose}. For
##' \code{windRose} and the \code{ws.int} default of 2 m/s, the
##' default, 4, generates the break points 2, 4, 6, 8 m/s. For
##' \code{pollutionRose}, the default, 6, attempts to breaks the
##' supplied data at approximately 6 sensible break points. However,
##' \code{breaks} can also be used to set specific break points. For
##' example, the argument \code{breaks = c(0, 1, 10, 100)} breaks the
##' data into segments <1, 1-10, 10-100, >100.
##' @param offset The size of the 'hole' in the middle of the plot, expressed
##' as a percentage of the polar axis scale, default 10.
##' @param normalise If \code{TRUE} each wind direction segment of a
##' pollution rose is normalised to equal one. This is useful for
##' showing how the concentrations (or other parameters) contribute to
##' each wind sector when the proprtion of time the wind is from that
##' direction is low. A line showing the probability that the wind
##' directions is from a particular wind sector is also shown.
##' @param max.freq Controls the scaling used by setting the maximum
##' value for the radial limits. This is useful to ensure several
##' plots use the same radial limits.
##' @param paddle Either \code{TRUE} (default) or \code{FALSE}. If \code{TRUE}
##' plots rose using `paddle' style spokes. If \code{FALSE} plots rose using
##' `wedge' style spokes.
##' @param key.header Adds additional text/labels above and/or below
##' the scale key, respectively. For example, passing
##' \code{windRose(mydata, key.header = "ws")} adds the addition text
##' as a scale header. Note: This argument is passed to
##' \code{drawOpenKey} via \code{quickText}, applying the auto.text
##' argument, to handle formatting.
##' @param key.footer see \code{key.footer}.
##' @param key.position Location where the scale key is to plotted.
##' Allowed arguments currently include \dQuote{top},
##' \dQuote{right}, \dQuote{bottom} and \dQuote{left}.
##' @param key Fine control of the scale key via \code{drawOpenKey}. See
##' \code{drawOpenKey} for further details.
##' @param dig.lab The number of signficant figures at which scientific number
##' formatting is used in break point and key labelling. Default 5.
##' @param statistic The \code{statistic} to be applied to each data
##' bin in the plot. Options currently include \dQuote{prop.count},
##' \dQuote{prop.mean} and \dQuote{abs.count}. The default
##' \dQuote{prop.count} sizes bins according to the proportion of
##' the frequency of measurements. Similarly, \dQuote{prop.mean} sizes
##' bins according to their relative contribution to the
##' mean. \dQuote{abs.count} provides the absolute count of
##' measurements in each bin.
##' @param pollutant Alternative data series to be sampled instead of wind
##' speed. The \code{windRose} default NULL is equivalent to \code{pollutant
##' = "ws"}.
##' @param annotate If \code{TRUE} then the percentage calm and mean
##' values are printed in each panel together with a description of
##' the statistic below the plot.
##' @param angle.scale The wind speed scale is by default shown at a
##' 315 degree angle. Sometimes the placement of the scale may
##' interfere with an interesting feature. The user can therefore set
##' \code{angle.scale} to another value (between 0 and 360 degrees) to
##' mitigate such problems. For example \code{angle.scale = 45} will
##' draw the scale heading in a NE direction.
##' @param border Border colour for shaded areas. Default is no border.
##' @param ... For \code{pollutionRose} other parameters that are
##' passed on to \code{windRose}. For \code{windRose} other parameters
##' that are passed on to \code{drawOpenKey}, \code{lattice:xyplot}
##' and \code{cutData}. Axis and title labelling options (\code{xlab},
##' \code{ylab}, \code{main}) are passed to \code{xyplot} via
##' \code{quickText} to handle routine formatting.
##'
##' @export windRose pollutionRose
##' @useDynLib ggopenair
##' @import dplyr
##' @import lazyeval
##' @importFrom plyr ddply ldply dlply llply numcolwise .
##' @importFrom graphics abline
##' @importFrom grDevices col2rgb colorRampPalette grey rgb xy.coords
##' @importFrom methods is
##' @importFrom stats aggregate approx as.dendrogram as.dist
##' ave coef cor dist formula hclust lm median
##' na.omit optimize order.dendrogram predict qchisq qnorm
##' qt quantile reshape sd smooth.spline spline stl ts
##' update var
##' @importFrom utils compareVersion modifyList packageDescription
##' read.csv read.table
##' @return As well as generating the plot itself, \code{windRose} and
##' \code{pollutionRose} also return an object of class \dQuote{openair}. The
##' object includes three main components: \code{call}, the command used to
##' generate the plot; \code{data}, the data frame of summarised information
##' used to make the plot; and \code{plot}, the plot itself. If retained, e.g.
##' using \code{output <- windRose(mydata)}, this output can be used to
##' recover the data, reproduce or rework the original plot or
##' undertake further analysis.
##'
##' An openair output can be manipulated using a number of generic operations,
##' including \code{print}, \code{plot} and \code{summarise}.
##'
##' Summarised proportions can also be extracted directly using the
##' \code{$data} operator, e.g. \code{object$data} for \code{output <-
##' windRose(mydata)}. This returns a data frame with three set columns:
##' \code{cond}, conditioning based on \code{type}; \code{wd}, the wind
##' direction; and \code{calm}, the \code{statistic} for the proportion of
##' data unattributed to any specific wind direction because it was collected
##' under calm conditions; and then several (one for each range binned for
##' the plot) columns giving proportions of measurements associated with each
##' \code{ws} or \code{pollutant} range plotted as a discrete panel.
##' @note \code{windRose} and \code{pollutionRose} both use \link{drawOpenKey}
##' to produce scale keys.
##' @author David Carslaw (with some additional contributions by Karl Ropkins)
##' @seealso See \code{\link{drawOpenKey}} for fine control of the scale key.
##'
##' See \code{\link{polarFreq}} for a more flexible version that considers
##' other statistics and pollutant concentrations.
##' @keywords methods
##' @references
##'
##' Applequist, S, 2012: Wind Rose Bias
##' Correction. J. Appl. Meteor. Climatol., 51, 1305-1309.
##'
##' This paper seems to be the original?
##'
##' Droppo, J.G. and B.A. Napier (2008) Wind Direction Bias in
##' Generating Wind Roses and Conducting Sector-Based Air Dispersion
##' Modeling, Journal of the Air & Waste Management Association, 58:7, 913-918.
##'
##' @examples
##'
##' # load example data from package data(mydata)
##'
##' # basic plot
##' windRose(mydata)
##'
##' # one windRose for each year
##' windRose(mydata,type = "year")
##'
##' # windRose in 10 degree intervals with gridlines and width adjusted
##' \dontrun{
##' windRose(mydata, angle = 10, width = 0.2, grid.line = 1)
##' }
##'
##' # pollutionRose of nox
##' pollutionRose(mydata, pollutant = "nox")
##'
##' ## source apportionment plot - contribution to mean
##' \dontrun{
##' pollutionRose(mydata, pollutant = "pm10", type = "year", statistic = "prop.mean")
##' }
##'
##' ## example of comparing 2 met sites
##' ## first we will make some new ws/wd data with a postive bias
##' mydata$ws2 = mydata$ws + 2 * rnorm(nrow(mydata)) + 1
##' mydata$wd2 = mydata$wd + 30 * rnorm(nrow(mydata)) + 30
##'
##' ## need to correct negative wd
##' id <- which(mydata$wd2 < 0)
##' mydata$wd2[id] <- mydata$wd2[id] + 360
##'
##' ## results show postive bias in wd and ws
##' pollutionRose(mydata, ws = "ws", wd = "wd", ws2 = "ws2", wd2 = "wd2")
windRose <- function (mydata, ws = "ws", wd = "wd", ws2 = NA, wd2 = NA,
ws.int = 2, angle = 30, type = "default", bias.corr = TRUE,
cols = "default", grid.line = NULL, width = 1, seg = NULL,
auto.text = TRUE, breaks = 4, offset = 10, normalise = FALSE,
max.freq = NULL, paddle = TRUE, key.header = NULL,
key.footer = "(m/s)", key.position = "bottom",
key = TRUE, dig.lab = 5, statistic = "prop.count",
pollutant = NULL, annotate = TRUE, angle.scale = 315, border = NA,
...)
{
if (is.null(seg)) seg <- 0.9
# make sure ws and wd and numeric
mydata <- checkNum(mydata, vars = c(ws, wd))
if (360 / angle != round(360 / angle)) {
warning("In windRose(...):\n angle will produce some spoke overlap",
"\n suggest one of: 5, 6, 8, 9, 10, 12, 15, 30, 45, etc.", call. = FALSE)
}
if (angle < 3) {
warning("In windRose(...):\n angle too small",
"\n enforcing 'angle = 3'", call. = FALSE)
angle <- 3
}
## extra args setup
extra <- list(...)
## label controls
extra$xlab <- if("xlab" %in% names(extra))
quickText(extra$xlab, auto.text) else quickText("", auto.text)
extra$ylab <- if("ylab" %in% names(extra))
quickText(extra$ylab, auto.text) else quickText("", auto.text)
extra$main <- if("main" %in% names(extra))
quickText(extra$main, auto.text) else quickText("", auto.text)
rounded <- FALSE ## is the wd already rounded to 10 degrees, if so need to correct bias later
if (all(mydata[[wd]] %% 10 == 0, na.rm = TRUE)) rounded <- TRUE
## preset statitistics
if (is.character(statistic)) {
## allowed cases
ok.stat <- c("prop.count", "prop.mean", "abs.count", "frequency")
if (!is.character(statistic) || !statistic[1] %in% ok.stat) {
warning("In windRose(...):\n statistic unrecognised",
"\n enforcing statistic = 'prop.count'", call. = FALSE)
statistic <- "prop.count"
}
if (statistic == "prop.count"){
stat.fun <- length
stat.unit <- "%"
stat.scale <- "all"
stat.lab <- "Frequency of counts by wind direction (%)"
stat.fun2 <- function(x) format(mean(x, na.rm = TRUE), digits = 5)
stat.lab2 <- "mean"
stat.labcalm <- function(x) round(x, 1)
}
if (statistic == "prop.mean") {
stat.fun <- function(x) sum(x, na.rm = TRUE)
stat.unit <- "%"
stat.scale <- "panel"
stat.lab <- "Proportion contribution to the mean (%)"
stat.fun2 <- function(x) format(mean(x, na.rm = TRUE), digits = 5)
stat.lab2 <- "mean"
stat.labcalm <- function(x) round(x, 1)
}
if (statistic == "abs.count" | statistic == "frequency") {
stat.fun <- length
stat.unit <- ""
stat.scale <- "none"
stat.lab <- "Count by wind direction"
stat.fun2 <- function(x) round(length(x), 0)
stat.lab2 <- "count"
stat.labcalm <- function(x) round(x, 0)
}
}
## variables we need
vars <- c(wd, ws)
diff <- FALSE ## i.e. not two sets of ws/wd
rm.neg <- TRUE ## will remove negative ws in check.prep
## case where two met data sets are to be compared
if (!is.na(ws2) & !is.na(wd2)) {
vars <- c(vars, ws2, wd2)
diff <- TRUE
rm.neg <- FALSE
mydata$ws <- mydata[[ws2]] - mydata[[ws]]
mydata$wd <- mydata[[wd2]] - mydata[[wd]]
## fix negative wd
id <- which(mydata$wd < 0)
if (length(id) > 0) mydata$wd[id] <- mydata$wd[id] + 360
pollutant <- "ws"
key.footer <- "ws"
wd <- "wd" ; ws <- "ws"
vars <- c("ws", "wd")
if (missing(angle)) angle <- 10
if (missing(offset)) offset <- 20
## set the breaks to cover all the data
if (is.na(breaks[1])) {
max.br <- max(ceiling(abs(c(min(mydata$ws, na.rm = TRUE),
max(mydata$ws, na.rm = TRUE)))))
breaks <- c(-1 * max.br, 0, max.br)
}
if (missing(cols)) cols <- c("lightskyblue", "tomato")
seg <- 1
}
if (any(type %in% dateTypes)) vars <- c(vars, "date")
if (!is.null(pollutant)) vars <- c(vars, pollutant)
mydata <- cutData(mydata, type, ...)
mydata <- checkPrep(mydata, vars, type, remove.calm = FALSE, remove.neg = rm.neg)
# can remove all missing data
mydata <- na.omit(mydata)
# remove lines where ws is missing
# wd can be NA and ws 0 (calm)
id <- which(is.na(mydata[[ws]]))
if (length(id) > 0)
mydata <- mydata[-id, ]
if (is.null(pollutant)) pollutant <- ws
mydata$x <- mydata[[pollutant]]
mydata[[wd]] <- angle * ceiling(mydata[[wd]] / angle - 0.5)
mydata[[wd]][mydata[[wd]] == 0] <- 360
## flag calms as negatives
mydata[[wd]][mydata[ , ws] == 0] <- -999 ## set wd to flag where there are calms
## do after rounding or -999 changes
if (length(breaks) == 1) breaks <- 0:(breaks - 1) * ws.int
if (max(breaks) < max(mydata$x, na.rm = TRUE))
breaks <- c(breaks, max(mydata$x, na.rm = TRUE))
if (min(breaks) > min(mydata$x, na.rm = TRUE))
warning ("Some values are below minimum break.")
breaks <- unique(breaks)
mydata$x <- cut(mydata$x, breaks = breaks, include.lowest = FALSE,
dig.lab = dig.lab)
## clean up cut intervals
labs <- gsub("[(]|[)]|[[]|[]]", "", levels(mydata$x))
labs <- gsub("[,]", " to ", labs)
## statistic handling
prepare.grid <- function(mydata) {
## these are all calms...
if (all(is.na(mydata$x))) {
weights <- data.frame(Interval1 = NA, wd = NA,
calm = 100, panel.fun = NA, mean.wd = NA, freqs = NA)
} else {
levels(mydata$x) <- c(paste("Interval", 1:length(labs), sep = ""))
all <- stat.fun(mydata[[wd]])
calm <- mydata[mydata[[wd]] == -999, ][[pollutant]]
calm <- stat.fun(calm)
weights <- tapply(mydata[[pollutant]], list(mydata[[wd]], mydata$x),
stat.fun)
freqs <- tapply(mydata[[pollutant]], mydata[[wd]], length)
## scaling
if (stat.scale == "all") {
calm <- calm / all
weights <- weights / all
}
if (stat.scale == "panel") {
temp <- stat.fun(stat.fun(weights)) + calm
calm <- calm / temp
weights <- weights / temp
}
weights[is.na(weights)] <- 0
weights <- t(apply(weights, 1, cumsum))
if (stat.scale == "all" | stat.scale == "panel"){
weights <- weights * 100
calm <- calm * 100
}
panel.fun <- stat.fun2(mydata[[pollutant]])
## calculate mean wd - useful for cases comparing two met data sets
u <- mean(sin(2 * pi * mydata[[wd]] / 360))
v <- mean(cos(2 * pi * mydata[[wd]] / 360))
mean.wd <- atan2(u, v) * 360 / 2 / pi
if (all(is.na(mean.wd))) {
mean.wd <- NA
} else {
if (mean.wd < 0) mean.wd <- mean.wd + 360
## show as a negative (bias)
if (mean.wd > 180) mean.wd <- mean.wd - 360
}
weights <- bind_cols(as_data_frame(weights),
data_frame(wd = as.numeric(row.names(weights)),
calm = calm, panel.fun = panel.fun,
mean.wd = mean.wd, freqs = freqs))
}
weights
}
if (paddle) {
poly <- function(wd, len1, len2, width, colour, x.off = 0, y.off = 0) {
theta <- wd * pi / 180
len1 <- len1 + off.set
len2 <- len2 + off.set
x1 <- len1 * sin(theta) - width * cos(theta) + x.off
x2 <- len1 * sin(theta) + width * cos(theta) + x.off
x3 <- len2 * sin(theta) - width * cos(theta) + x.off
x4 <- len2 * sin(theta) + width * cos(theta) + x.off
y1 <- len1 * cos(theta) + width * sin(theta) + y.off
y2 <- len1 * cos(theta) - width * sin(theta) + y.off
y3 <- len2 * cos(theta) + width * sin(theta) + y.off
y4 <- len2 * cos(theta) - width * sin(theta) + y.off
lpolygon(c(x1, x2, x4, x3), c(y1, y2, y4, y3), col = colour,
border = border)
}
} else {
poly <- function(wd, len1, len2, width, colour, x.off = 0,
y.off = 0) {
len1 <- len1 + off.set
len2 <- len2 + off.set
theta <- seq((wd - seg * angle / 2), (wd + seg * angle / 2),
length.out = (angle - 2) * 10)
theta <- ifelse(theta < 1, 360 - theta, theta)
theta <- theta * pi / 180
x1 <- len1 * sin(theta) + x.off
x2 <- rev(len2 * sin(theta) + x.off)
y1 <- len1 * cos(theta) + x.off
y2 <- rev(len2 * cos(theta) + x.off)
lpolygon(c(x1, x2), c(y1, y2), col = colour, border = border)
}
}
results <- group_by_(mydata, .dots = type) %>%
do(prepare.grid(.))
## format
results$calm <- stat.labcalm(results$calm)
results$mean.wd <- stat.labcalm(results$mean.wd)
## correction for bias when angle does not divide exactly into 360
if (bias.corr & rounded) {
wd <- seq(10, 360, 10)
tmp <- angle * ceiling(wd / angle - 0.5)
id <- which(tmp == 0)
if (length(id > 0)) tmp[id] <- 360
tmp <- table(tmp) ## number of sectors spanned
vars <- grep("Interval[1-9]", names(results)) ## the frequencies
results[-1, vars] <- results[-1, vars] * mean(tmp) /tmp
}
## proper names of labelling###########################################
strip.dat <- strip.fun(results, type, auto.text)
pol.name <- strip.dat[[1]]
pol.name2 <- strip.dat[[2]]
# proper names for strip labels
levels(results[[type[1]]]) <- pol.name
if (length(type) == 2L) levels(results[[type[2]]]) <- pol.name2
if (length(labs) < length(cols)) {
col <- cols[1:length(labs)]
} else {
col <- openColours(cols, length(labs))
}
## normalise by sector
if (normalise) {
vars <- grep("Interval[1-9]", names(results))
## original frequencies, so we can plot the wind frequency line
results$freq <- results[[max(vars)]]
results$freq <- ave(results$freq, results[type], FUN = function(x) x / sum(x))
## scale by maximum frequency
results$norm <- results$freq / max(results$freq)
## normalise
results[[vars]] <- results[[vars]] / results[[max(vars)]]
stat.lab <- "Normalised by wind sector"
stat.unit <- ""
}
if (is.null(max.freq)) {
max.freq <- max(results[, (length(type) + 1):(length(labs) +
length(type))],
na.rm = TRUE)
} else {
max.freq <- max.freq
}
off.set <- max.freq * (offset / 100)
box.widths <- seq(0.002 ^ 0.25, 0.016 ^ 0.25,
length.out = length(labs)) ^ 4
box.widths <- box.widths * max.freq * angle / 5
mymax <- 2 * max.freq
myby <- if (is.null(grid.line)) pretty(c(0, mymax), 10)[2] else grid.line
if (myby / mymax > 0.9) myby <- mymax * 0.9
if (annotate) sub <- stat.lab else sub <- NULL
# make a data frame of widths, n is the number of the interval
intervals <- names(results[grep(pattern = "Interval", names(results))])
intervals <- data_frame(interval = intervals, width = box.widths,
n = seq_along(intervals))
# add Interval0
results$Interval0 <- 0
# reshape for calculations
results <- gather(results, key = interval, value = value, contains("Interval"))
# only need some columns for this
results <- select_(results, .dots = c(type, "wd", "interval", "value"))
results <- filter(results, wd != -999) # don't need calms to plot polygons
# add widths
results <- full_join(results, intervals, by = "interval")
results <- mutate(results, n = ifelse(is.na(n), 0, n))
# these results are used for plotting
results2 <- group_by_(results, .dots = c(type, "wd")) %>%
do(calc_poly_coord(., off.set))
angles <- seq(0, 2 * pi, length = 360)
# data for underlying background
datapoly <- data.frame(x = (max.freq + off.set) * sin(angles),
y = (max.freq + off.set) * cos(angles))
# data for gridlines
lines_at <- seq(off.set, mymax, by = myby)
datalines <- data.frame(value = rep(lines_at, each = 360),
x = rep(lines_at, each = 360) * sin(angles),
y = rep(lines_at, each = 360) * cos(angles))
# data for axes
upper <- max.freq + off.set
data.axes <- data.frame(x1 = c(-upper, 0), y1 = c(0, -upper),
x2 = c(upper, 0), y2 = c(0, upper))
plt <- ggplot(results2, aes(x, y, group = interaction(interval, wd))) +
theme_void() +
coord_cartesian(xlim = 1.03 * c(-max.freq - off.set, max.freq + off.set),
ylim = 1.03 * c(-max.freq - off.set, max.freq + off.set)) +
geom_polygon(data = datapoly, aes(x = x, y = y), color = "grey95",
fill = "grey95", inherit.aes = FALSE) +
geom_path(data = datalines, aes(x = x, y = y), color = "white",
inherit.aes = FALSE) +
geom_segment(data = data.axes, aes(x = x1, y = y1, xend = x2, yend = y2),
color = "white", inherit.aes = FALSE) +
annotate("text", upper * -1 * 0.95, 0, label = "W") +
annotate("text", 0, upper * -1 * 0.95, label = "S") +
annotate("text", 0, upper * 0.95, label = "N") +
annotate("text", upper * 0.95, 0, label = "E") +
geom_polygon(aes(fill = interval)) +
scale_fill_manual(values = openColours(cols, max(results$n)),
labels = labs) +
theme(aspect.ratio = 1, strip.background = element_rect(fill = "grey"))
if (!"default" %in% type) {
if (length(type) == 1L)
plt <- plt + facet_wrap(reformulate(type), labeller = label_parsed)
if (length(type) == 2L)
plt <- plt + facet_grid(paste(type[2], "~", type[1]), labeller = label_parsed)
}
plot(plt)
newdata <- results
output <- list(plot = plt, data = newdata, call = match.call())
class(output) <- "openair"
invisible(output)
}
## adds a line showing probability wind direction is from a particular sector
## used when normalise = TRUE
panel.wdprob <- function(dat, seg, angle, off.set) {
len1 <- off.set
x.off <- 0; y.off <- 0
makeline <- function(i, dat) {
theta <- seq((dat$wd[i] - seg * angle / 2), (dat$wd[i] + seg * angle / 2),
length.out = (angle - 2) * 10)
theta <- ifelse(theta < 1, 360 - theta, theta)
theta <- theta * pi / 180
x1 <- len1 * sin(theta) + x.off
x2 <- rev((dat$norm[i] + off.set) * sin(theta) + x.off)
y1 <- len1 * cos(theta) + x.off
y2 <- rev((dat$norm[i] + off.set) * cos(theta) + x.off)
lpolygon(c(x1, x2), c(y1, y2), col = "transparent", border = "black", lwd = 2)
}
lapply(1:nrow(dat), makeline, dat)
}
calc_poly_coord <- function(data, off.set) {
# order by n
data <- arrange(data, n)
data$wd <- data$wd * pi / 180
data$value <- data$value + off.set
# groups of 4 points to make a polygon, then back to starting point
coord_x <- sapply(2:(nrow(data)),
function (i)
with(data,
c(value[i - 1] * sin(wd[i - 1]) -
width[i] * cos(wd[i - 1]),
value[i - 1] * sin(wd[i - 1]) +
width[i] * cos(wd[i - 1]),
value[i] * sin(wd[i]) +
width[i] * cos(wd[i]),
value[i] * sin(wd[i]) -
width[i] * cos(wd[i]),
value[i - 1] * sin(wd[i - 1]) -
width[i] * cos(wd[i - 1]))))
coord_x <- as.vector(coord_x)
coord_y <- sapply(2:(nrow(data)),
function (i)
with(data,
c(value[i - 1] * cos(wd[i - 1]) +
width[i] * sin(wd[i - 1]),
value[i - 1] * cos(wd[i - 1]) -
width[i] * sin(wd[i - 1]),
value[i] * cos(wd[i]) -
width[i] * sin(wd[i]),
value[i] * cos(wd[i]) +
width[i] * sin(wd[i]),
value[i - 1] * cos(wd[i - 1]) +
width[i] * sin(wd[i - 1]))))
coord_y <- as.vector(coord_y)
# make a data frame with all the coordinates for the polygons
result <- data_frame(x = coord_x,
y = coord_y,
interval = rep(data$interval[2:nrow(data)],
each = 5))
return(result)
}
davidcarslaw/ggopenair documentation built on May 14, 2019, 10:37 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
pollutionRose <- function(mydata, pollutant = "nox", key.footer = pollutant,
key.position = "right", key = TRUE,
breaks = 6, paddle = FALSE, seg = 0.9, normalise = FALSE,
...)
{
## extra args setup
extra <- list(...)
## check to see if two met data sets are being compared.
## if so, set pollutant to one of the names
if ("ws2" %in% names(extra)) {
pollutant <- extra$ws
if (missing(breaks)) breaks <- NA
}
if (is.null(breaks)) breaks <- 6
if (is.numeric(breaks) & length(breaks) == 1) {
## breaks from the minimum to 90th percentile, which generally gives sensible
## spacing for skewed data. Maximum is added later.
breaks <- unique(pretty(c(min(mydata[[pollutant]], na.rm = TRUE),
quantile(mydata[[pollutant]], probs = 0.9, na.rm = TRUE),
breaks)))
}
windRose(mydata, pollutant = pollutant, paddle = paddle, seg = seg,
key.position = key.position, key.footer = key.footer, key = key,
breaks = breaks, normalise = normalise, ...)
}
##' Traditional wind rose plot and pollution rose variation
##'
##' The traditional wind rose plot that plots wind speed and wind direction by
##' different intervals. The pollution rose applies the same plot structure but
##' substitutes other measurements, most commonly a pollutant time series, for
##' wind speed.
##'
##' For \code{windRose} data are summarised by direction, typically by 45 or 30
##' (or 10) degrees and by different wind speed categories. Typically, wind
##' speeds are represented by different width "paddles". The plots show the
##' proportion (here represented as a percentage) of time that the wind is from
##' a certain angle and wind speed range.
##'
##' By default \code{windRose} will plot a windRose in using "paddle" style
##' segments and placing the scale key below the plot.
##'
##' The argument \code{pollutant} uses the same plotting structure but
##' substitutes another data series, defined by \code{pollutant}, for wind
##' speed.
##'
##' The option \code{statistic = "prop.mean"} provides a measure of the
##' relative contribution of each bin to the panel mean, and is intended for
##' use with \code{pollutionRose}.
##'
##' \code{pollutionRose} is a \code{windRose} wrapper which brings
##' \code{pollutant} forward in the argument list, and attempts to sensibly
##' rescale break points based on the \code{pollutant} data range by by-passing
##' \code{ws.int}.
##'
##' By default, \code{pollutionRose} will plot a pollution rose of \code{nox}
##' using "wedge" style segments and placing the scale key to the right of the
##' plot.
##'
##' It is possible to compare two wind speed-direction data sets using
##' \code{pollutionRose}. There are many reasons for doing so e.g. to
##' see how one site compares with another or for meteorological model
##' evaluation. In this case, \code{ws} and \code{wd} are considered
##' to the the reference data sets with which a second set of wind
##' speed and wind directions are to be compared (\code{ws2} and
##' \code{wd2}). The first set of values is subtracted from the second
##' and the differences compared. If for example, \code{wd2} was
##' biased positive compared with \code{wd} then \code{pollutionRose}
##' will show the bias in polar coordinates. In its default use, wind
##' direction bias is colour-coded to show negative bias in one colour
##' and positive bias in another.
##'
##' @usage windRose(mydata, ws = "ws", wd = "wd", ws2 = NA, wd2 = NA,
##' ws.int = 2, angle = 30, type = "default", bias.corr = TRUE, cols = "default",
##' grid.line = NULL, width = 1, seg = NULL, auto.text = TRUE, breaks
##' = 4, offset = 10, normalise = FALSE, max.freq = NULL, paddle = TRUE, key.header =
##' NULL, key.footer = "(m/s)", key.position = "bottom", key = TRUE,
##' dig.lab = 5, statistic = "prop.count", pollutant = NULL, annotate
##' = TRUE, angle.scale = 315, border = NA, ...)
##'
##'
##' pollutionRose(mydata, pollutant = "nox", key.footer = pollutant,
##' key.position = "right", key = TRUE, breaks = 6, paddle = FALSE,
##' seg = 0.9, normalise = FALSE, ...)
##'
##'
##' @aliases windRose pollutionRose
##' @param mydata A data frame containing fields \code{ws} and \code{wd}
##' @param ws Name of the column representing wind speed.
##' @param wd Name of the column representing wind direction.
##' @param ws2 The user can supply a second set of wind speed and wind
##' direction values with which the first can be compared. See details
##' below for full explanation.
##' @param wd2 see \code{ws2}.
##' @param ws.int The Wind speed interval. Default is 2 m/s but for low met
##' masts with low mean wind speeds a value of 1 or 0.5 m/s may be better.
##' Note, this argument is superseded in \code{pollutionRose}. See
##' \code{breaks} below.
##' @param angle Default angle of \dQuote{spokes} is 30. Other potentially useful
##' angles are 45 and 10. Note that the width of the wind speed interval may
##' need adjusting using \code{width}.
##' @param type \code{type} determines how the data are split
##' i.e. conditioned, and then plotted. The default is will produce a
##' single plot using the entire data. Type can be one of the built-in
##' types as detailed in \code{cutData} e.g. \dQuote{season},
##' \dQuote{year}, \dQuote{weekday} and so on. For example, \code{type
##' = "season"} will produce four plots --- one for each season.
##'
##' It is also possible to choose \code{type} as another variable in the data
##' frame. If that variable is numeric, then the data will be split into four
##' quantiles (if possible) and labelled accordingly. If type is an existing
##' character or factor variable, then those categories/levels will be used
##' directly. This offers great flexibility for understanding the variation
##' of different variables and how they depend on one another.
##'
##' Type can be up length two e.g. \code{type = c("season", "weekday")} will
##' produce a 2x2 plot split by season and day of the week. Note, when two
##' types are provided the first forms the columns and the second the rows.
##' @param bias.corr When \code{angle} does not divide exactly into
##' 360 a bias is introduced in the frequencies when the wind
##' direction is already supplied rounded to the nearest 10 degrees,
##' as is often the case. For example, if \code{angle = 22.5}, N, E,
##' S, W will include 3 wind sectors and all other angles will be
##' two. A bias correction can made to correct for this problem. A
##' simple method according to Applequist (2012) is used to adjust the
##' frequencies.
##' @param cols Colours to be used for plotting. Options include
##' \dQuote{default}, \dQuote{increment}, \dQuote{heat}, \dQuote{jet},
##' \dQuote{hue} and user defined. For user defined the user can
##' supply a list of colour names recognised by R (type
##' \code{colours()} to see the full list). An example would be
##' \code{cols = c("yellow", "green", "blue", "black")}.
##' @param grid.line Grid line interval to use. If \code{NULL}, as in default,
##' this is assigned by \code{windRose} based on the available data range.
##' However, it can also be forced to a specific value, e.g.
##' \code{grid.line = 10}.
##' @param width For \code{paddle = TRUE}, the adjustment factor for width of
##' wind speed intervals. For example, \code{width = 1.5} will make the
##' paddle width 1.5 times wider.
##' @param seg For \code{pollutionRose} \code{seg} determines with
##' width of the segments. For example, \code{seg = 0.5} will produce
##' segments 0.5 * \code{angle}.
##' @param auto.text Either \code{TRUE} (default) or \code{FALSE}. If
##' \code{TRUE} titles and axis labels will automatically try and
##' format pollutant names and units properly e.g. by subscripting
##' the \sQuote{2} in NO2.
##' @param breaks Most commonly, the number of break points for wind
##' speed in \code{windRose} or pollutant in \code{pollutionRose}. For
##' \code{windRose} and the \code{ws.int} default of 2 m/s, the
##' default, 4, generates the break points 2, 4, 6, 8 m/s. For
##' \code{pollutionRose}, the default, 6, attempts to breaks the
##' supplied data at approximately 6 sensible break points. However,
##' \code{breaks} can also be used to set specific break points. For
##' example, the argument \code{breaks = c(0, 1, 10, 100)} breaks the
##' data into segments <1, 1-10, 10-100, >100.
##' @param offset The size of the 'hole' in the middle of the plot, expressed
##' as a percentage of the polar axis scale, default 10.
##' @param normalise If \code{TRUE} each wind direction segment of a
##' pollution rose is normalised to equal one. This is useful for
##' showing how the concentrations (or other parameters) contribute to
##' each wind sector when the proprtion of time the wind is from that
##' direction is low. A line showing the probability that the wind
##' directions is from a particular wind sector is also shown.
##' @param max.freq Controls the scaling used by setting the maximum
##' value for the radial limits. This is useful to ensure several
##' plots use the same radial limits.
##' @param paddle Either \code{TRUE} (default) or \code{FALSE}. If \code{TRUE}
##' plots rose using `paddle' style spokes. If \code{FALSE} plots rose using
##' `wedge' style spokes.
##' @param key.header Adds additional text/labels above and/or below
##' the scale key, respectively. For example, passing
##' \code{windRose(mydata, key.header = "ws")} adds the addition text
##' as a scale header. Note: This argument is passed to
##' \code{drawOpenKey} via \code{quickText}, applying the auto.text
##' argument, to handle formatting.
##' @param key.footer see \code{key.footer}.
##' @param key.position Location where the scale key is to plotted.
##' Allowed arguments currently include \dQuote{top},
##' \dQuote{right}, \dQuote{bottom} and \dQuote{left}.
##' @param key Fine control of the scale key via \code{drawOpenKey}. See
##' \code{drawOpenKey} for further details.
##' @param dig.lab The number of signficant figures at which scientific number
##' formatting is used in break point and key labelling. Default 5.
##' @param statistic The \code{statistic} to be applied to each data
##' bin in the plot. Options currently include \dQuote{prop.count},
##' \dQuote{prop.mean} and \dQuote{abs.count}. The default
##' \dQuote{prop.count} sizes bins according to the proportion of
##' the frequency of measurements. Similarly, \dQuote{prop.mean} sizes
##' bins according to their relative contribution to the
##' mean. \dQuote{abs.count} provides the absolute count of
##' measurements in each bin.
##' @param pollutant Alternative data series to be sampled instead of wind
##' speed. The \code{windRose} default NULL is equivalent to \code{pollutant
##' = "ws"}.
##' @param annotate If \code{TRUE} then the percentage calm and mean
##' values are printed in each panel together with a description of
##' the statistic below the plot.
##' @param angle.scale The wind speed scale is by default shown at a
##' 315 degree angle. Sometimes the placement of the scale may
##' interfere with an interesting feature. The user can therefore set
##' \code{angle.scale} to another value (between 0 and 360 degrees) to
##' mitigate such problems. For example \code{angle.scale = 45} will
##' draw the scale heading in a NE direction.
##' @param border Border colour for shaded areas. Default is no border.
##' @param ... For \code{pollutionRose} other parameters that are
##' passed on to \code{windRose}. For \code{windRose} other parameters
##' that are passed on to \code{drawOpenKey}, \code{lattice:xyplot}
##' and \code{cutData}. Axis and title labelling options (\code{xlab},
##' \code{ylab}, \code{main}) are passed to \code{xyplot} via
##' \code{quickText} to handle routine formatting.
##'
##' @export windRose pollutionRose
##' @useDynLib ggopenair
##' @import dplyr
##' @import lazyeval
##' @importFrom plyr ddply ldply dlply llply numcolwise .
##' @importFrom graphics abline
##' @importFrom grDevices col2rgb colorRampPalette grey rgb xy.coords
##' @importFrom methods is
##' @importFrom stats aggregate approx as.dendrogram as.dist
##' ave coef cor dist formula hclust lm median
##' na.omit optimize order.dendrogram predict qchisq qnorm
##' qt quantile reshape sd smooth.spline spline stl ts
##' update var
##' @importFrom utils compareVersion modifyList packageDescription
##' read.csv read.table
##' @return As well as generating the plot itself, \code{windRose} and
##' \code{pollutionRose} also return an object of class \dQuote{openair}. The
##' object includes three main components: \code{call}, the command used to
##' generate the plot; \code{data}, the data frame of summarised information
##' used to make the plot; and \code{plot}, the plot itself. If retained, e.g.
##' using \code{output <- windRose(mydata)}, this output can be used to
##' recover the data, reproduce or rework the original plot or
##' undertake further analysis.
##'
##' An openair output can be manipulated using a number of generic operations,
##' including \code{print}, \code{plot} and \code{summarise}.
##'
##' Summarised proportions can also be extracted directly using the
##' \code{$data} operator, e.g. \code{object$data} for \code{output <-
##' windRose(mydata)}. This returns a data frame with three set columns:
##' \code{cond}, conditioning based on \code{type}; \code{wd}, the wind
##' direction; and \code{calm}, the \code{statistic} for the proportion of
##' data unattributed to any specific wind direction because it was collected
##' under calm conditions; and then several (one for each range binned for
##' the plot) columns giving proportions of measurements associated with each
##' \code{ws} or \code{pollutant} range plotted as a discrete panel.
##' @note \code{windRose} and \code{pollutionRose} both use \link{drawOpenKey}
##' to produce scale keys.
##' @author David Carslaw (with some additional contributions by Karl Ropkins)
##' @seealso See \code{\link{drawOpenKey}} for fine control of the scale key.
##'
##' See \code{\link{polarFreq}} for a more flexible version that considers
##' other statistics and pollutant concentrations.
##' @keywords methods
##' @references
##'
##' Applequist, S, 2012: Wind Rose Bias
##' Correction. J. Appl. Meteor. Climatol., 51, 1305-1309.
##'
##' This paper seems to be the original?
##'
##' Droppo, J.G. and B.A. Napier (2008) Wind Direction Bias in
##' Generating Wind Roses and Conducting Sector-Based Air Dispersion
##' Modeling, Journal of the Air & Waste Management Association, 58:7, 913-918.
##'
##' @examples
##'
##' # load example data from package data(mydata)
##'
##' # basic plot
##' windRose(mydata)
##'
##' # one windRose for each year
##' windRose(mydata,type = "year")
##'
##' # windRose in 10 degree intervals with gridlines and width adjusted
##' \dontrun{
##' windRose(mydata, angle = 10, width = 0.2, grid.line = 1)
##' }
##'
##' # pollutionRose of nox
##' pollutionRose(mydata, pollutant = "nox")
##'
##' ## source apportionment plot - contribution to mean
##' \dontrun{
##' pollutionRose(mydata, pollutant = "pm10", type = "year", statistic = "prop.mean")
##' }
##'
##' ## example of comparing 2 met sites
##' ## first we will make some new ws/wd data with a postive bias
##' mydata$ws2 = mydata$ws + 2 * rnorm(nrow(mydata)) + 1
##' mydata$wd2 = mydata$wd + 30 * rnorm(nrow(mydata)) + 30
##'
##' ## need to correct negative wd
##' id <- which(mydata$wd2 < 0)
##' mydata$wd2[id] <- mydata$wd2[id] + 360
##'
##' ## results show postive bias in wd and ws
##' pollutionRose(mydata, ws = "ws", wd = "wd", ws2 = "ws2", wd2 = "wd2")
windRose <- function (mydata, ws = "ws", wd = "wd", ws2 = NA, wd2 = NA,
ws.int = 2, angle = 30, type = "default", bias.corr = TRUE,
cols = "default", grid.line = NULL, width = 1, seg = NULL,
auto.text = TRUE, breaks = 4, offset = 10, normalise = FALSE,
max.freq = NULL, paddle = TRUE, key.header = NULL,
key.footer = "(m/s)", key.position = "bottom",
key = TRUE, dig.lab = 5, statistic = "prop.count",
pollutant = NULL, annotate = TRUE, angle.scale = 315, border = NA,
...)
{
if (is.null(seg)) seg <- 0.9
# make sure ws and wd and numeric
mydata <- checkNum(mydata, vars = c(ws, wd))
if (360 / angle != round(360 / angle)) {
warning("In windRose(...):\n angle will produce some spoke overlap",
"\n suggest one of: 5, 6, 8, 9, 10, 12, 15, 30, 45, etc.", call. = FALSE)
}
if (angle < 3) {
warning("In windRose(...):\n angle too small",
"\n enforcing 'angle = 3'", call. = FALSE)
angle <- 3
}
## extra args setup
extra <- list(...)
## label controls
extra$xlab <- if("xlab" %in% names(extra))
quickText(extra$xlab, auto.text) else quickText("", auto.text)
extra$ylab <- if("ylab" %in% names(extra))
quickText(extra$ylab, auto.text) else quickText("", auto.text)
extra$main <- if("main" %in% names(extra))
quickText(extra$main, auto.text) else quickText("", auto.text)
rounded <- FALSE ## is the wd already rounded to 10 degrees, if so need to correct bias later
if (all(mydata[[wd]] %% 10 == 0, na.rm = TRUE)) rounded <- TRUE
## preset statitistics
if (is.character(statistic)) {
## allowed cases
ok.stat <- c("prop.count", "prop.mean", "abs.count", "frequency")
if (!is.character(statistic) || !statistic[1] %in% ok.stat) {
warning("In windRose(...):\n statistic unrecognised",
"\n enforcing statistic = 'prop.count'", call. = FALSE)
statistic <- "prop.count"
}
if (statistic == "prop.count"){
stat.fun <- length
stat.unit <- "%"
stat.scale <- "all"
stat.lab <- "Frequency of counts by wind direction (%)"
stat.fun2 <- function(x) format(mean(x, na.rm = TRUE), digits = 5)
stat.lab2 <- "mean"
stat.labcalm <- function(x) round(x, 1)
}
if (statistic == "prop.mean") {
stat.fun <- function(x) sum(x, na.rm = TRUE)
stat.unit <- "%"
stat.scale <- "panel"
stat.lab <- "Proportion contribution to the mean (%)"
stat.fun2 <- function(x) format(mean(x, na.rm = TRUE), digits = 5)
stat.lab2 <- "mean"
stat.labcalm <- function(x) round(x, 1)
}
if (statistic == "abs.count" | statistic == "frequency") {
stat.fun <- length
stat.unit <- ""
stat.scale <- "none"
stat.lab <- "Count by wind direction"
stat.fun2 <- function(x) round(length(x), 0)
stat.lab2 <- "count"
stat.labcalm <- function(x) round(x, 0)
}
}
## variables we need
vars <- c(wd, ws)
diff <- FALSE ## i.e. not two sets of ws/wd
rm.neg <- TRUE ## will remove negative ws in check.prep
## case where two met data sets are to be compared
if (!is.na(ws2) & !is.na(wd2)) {
vars <- c(vars, ws2, wd2)
diff <- TRUE
rm.neg <- FALSE
mydata$ws <- mydata[[ws2]] - mydata[[ws]]
mydata$wd <- mydata[[wd2]] - mydata[[wd]]
## fix negative wd
id <- which(mydata$wd < 0)
if (length(id) > 0) mydata$wd[id] <- mydata$wd[id] + 360
pollutant <- "ws"
key.footer <- "ws"
wd <- "wd" ; ws <- "ws"
vars <- c("ws", "wd")
if (missing(angle)) angle <- 10
if (missing(offset)) offset <- 20
## set the breaks to cover all the data
if (is.na(breaks[1])) {
max.br <- max(ceiling(abs(c(min(mydata$ws, na.rm = TRUE),
max(mydata$ws, na.rm = TRUE)))))
breaks <- c(-1 * max.br, 0, max.br)
}
if (missing(cols)) cols <- c("lightskyblue", "tomato")
seg <- 1
}
if (any(type %in% dateTypes)) vars <- c(vars, "date")
if (!is.null(pollutant)) vars <- c(vars, pollutant)
mydata <- cutData(mydata, type, ...)
mydata <- checkPrep(mydata, vars, type, remove.calm = FALSE, remove.neg = rm.neg)
# can remove all missing data
mydata <- na.omit(mydata)
# remove lines where ws is missing
# wd can be NA and ws 0 (calm)
id <- which(is.na(mydata[[ws]]))
if (length(id) > 0)
mydata <- mydata[-id, ]
if (is.null(pollutant)) pollutant <- ws
mydata$x <- mydata[[pollutant]]
mydata[[wd]] <- angle * ceiling(mydata[[wd]] / angle - 0.5)
mydata[[wd]][mydata[[wd]] == 0] <- 360
## flag calms as negatives
mydata[[wd]][mydata[ , ws] == 0] <- -999 ## set wd to flag where there are calms
## do after rounding or -999 changes
if (length(breaks) == 1) breaks <- 0:(breaks - 1) * ws.int
if (max(breaks) < max(mydata$x, na.rm = TRUE))
breaks <- c(breaks, max(mydata$x, na.rm = TRUE))
if (min(breaks) > min(mydata$x, na.rm = TRUE))
warning ("Some values are below minimum break.")
breaks <- unique(breaks)
mydata$x <- cut(mydata$x, breaks = breaks, include.lowest = FALSE,
dig.lab = dig.lab)
## clean up cut intervals
labs <- gsub("[(]|[)]|[[]|[]]", "", levels(mydata$x))
labs <- gsub("[,]", " to ", labs)
## statistic handling
prepare.grid <- function(mydata) {
## these are all calms...
if (all(is.na(mydata$x))) {
weights <- data.frame(Interval1 = NA, wd = NA,
calm = 100, panel.fun = NA, mean.wd = NA, freqs = NA)
} else {
levels(mydata$x) <- c(paste("Interval", 1:length(labs), sep = ""))
all <- stat.fun(mydata[[wd]])
calm <- mydata[mydata[[wd]] == -999, ][[pollutant]]
calm <- stat.fun(calm)
weights <- tapply(mydata[[pollutant]], list(mydata[[wd]], mydata$x),
stat.fun)
freqs <- tapply(mydata[[pollutant]], mydata[[wd]], length)
## scaling
if (stat.scale == "all") {
calm <- calm / all
weights <- weights / all
}
if (stat.scale == "panel") {
temp <- stat.fun(stat.fun(weights)) + calm
calm <- calm / temp
weights <- weights / temp
}
weights[is.na(weights)] <- 0
weights <- t(apply(weights, 1, cumsum))
if (stat.scale == "all" | stat.scale == "panel"){
weights <- weights * 100
calm <- calm * 100
}
panel.fun <- stat.fun2(mydata[[pollutant]])
## calculate mean wd - useful for cases comparing two met data sets
u <- mean(sin(2 * pi * mydata[[wd]] / 360))
v <- mean(cos(2 * pi * mydata[[wd]] / 360))
mean.wd <- atan2(u, v) * 360 / 2 / pi
if (all(is.na(mean.wd))) {
mean.wd <- NA
} else {
if (mean.wd < 0) mean.wd <- mean.wd + 360
## show as a negative (bias)
if (mean.wd > 180) mean.wd <- mean.wd - 360
}
weights <- bind_cols(as_data_frame(weights),
data_frame(wd = as.numeric(row.names(weights)),
calm = calm, panel.fun = panel.fun,
mean.wd = mean.wd, freqs = freqs))
}
weights
}
if (paddle) {
poly <- function(wd, len1, len2, width, colour, x.off = 0, y.off = 0) {
theta <- wd * pi / 180
len1 <- len1 + off.set
len2 <- len2 + off.set
x1 <- len1 * sin(theta) - width * cos(theta) + x.off
x2 <- len1 * sin(theta) + width * cos(theta) + x.off
x3 <- len2 * sin(theta) - width * cos(theta) + x.off
x4 <- len2 * sin(theta) + width * cos(theta) + x.off
y1 <- len1 * cos(theta) + width * sin(theta) + y.off
y2 <- len1 * cos(theta) - width * sin(theta) + y.off
y3 <- len2 * cos(theta) + width * sin(theta) + y.off
y4 <- len2 * cos(theta) - width * sin(theta) + y.off
lpolygon(c(x1, x2, x4, x3), c(y1, y2, y4, y3), col = colour,
border = border)
}
} else {
poly <- function(wd, len1, len2, width, colour, x.off = 0,
y.off = 0) {
len1 <- len1 + off.set
len2 <- len2 + off.set
theta <- seq((wd - seg * angle / 2), (wd + seg * angle / 2),
length.out = (angle - 2) * 10)
theta <- ifelse(theta < 1, 360 - theta, theta)
theta <- theta * pi / 180
x1 <- len1 * sin(theta) + x.off
x2 <- rev(len2 * sin(theta) + x.off)
y1 <- len1 * cos(theta) + x.off
y2 <- rev(len2 * cos(theta) + x.off)
lpolygon(c(x1, x2), c(y1, y2), col = colour, border = border)
}
}
results <- group_by_(mydata, .dots = type) %>%
do(prepare.grid(.))
## format
results$calm <- stat.labcalm(results$calm)
results$mean.wd <- stat.labcalm(results$mean.wd)
## correction for bias when angle does not divide exactly into 360
if (bias.corr & rounded) {
wd <- seq(10, 360, 10)
tmp <- angle * ceiling(wd / angle - 0.5)
id <- which(tmp == 0)
if (length(id > 0)) tmp[id] <- 360
tmp <- table(tmp) ## number of sectors spanned
vars <- grep("Interval[1-9]", names(results)) ## the frequencies
results[-1, vars] <- results[-1, vars] * mean(tmp) /tmp
}
## proper names of labelling###########################################
strip.dat <- strip.fun(results, type, auto.text)
pol.name <- strip.dat[[1]]
pol.name2 <- strip.dat[[2]]
# proper names for strip labels
levels(results[[type[1]]]) <- pol.name
if (length(type) == 2L) levels(results[[type[2]]]) <- pol.name2
if (length(labs) < length(cols)) {
col <- cols[1:length(labs)]
} else {
col <- openColours(cols, length(labs))
}
## normalise by sector
if (normalise) {
vars <- grep("Interval[1-9]", names(results))
## original frequencies, so we can plot the wind frequency line
results$freq <- results[[max(vars)]]
results$freq <- ave(results$freq, results[type], FUN = function(x) x / sum(x))
## scale by maximum frequency
results$norm <- results$freq / max(results$freq)
## normalise
results[[vars]] <- results[[vars]] / results[[max(vars)]]
stat.lab <- "Normalised by wind sector"
stat.unit <- ""
}
if (is.null(max.freq)) {
max.freq <- max(results[, (length(type) + 1):(length(labs) +
length(type))],
na.rm = TRUE)
} else {
max.freq <- max.freq
}
off.set <- max.freq * (offset / 100)
box.widths <- seq(0.002 ^ 0.25, 0.016 ^ 0.25,
length.out = length(labs)) ^ 4
box.widths <- box.widths * max.freq * angle / 5
mymax <- 2 * max.freq
myby <- if (is.null(grid.line)) pretty(c(0, mymax), 10)[2] else grid.line
if (myby / mymax > 0.9) myby <- mymax * 0.9
if (annotate) sub <- stat.lab else sub <- NULL
# make a data frame of widths, n is the number of the interval
intervals <- names(results[grep(pattern = "Interval", names(results))])
intervals <- data_frame(interval = intervals, width = box.widths,
n = seq_along(intervals))
# add Interval0
results$Interval0 <- 0
# reshape for calculations
results <- gather(results, key = interval, value = value, contains("Interval"))
# only need some columns for this
results <- select_(results, .dots = c(type, "wd", "interval", "value"))
results <- filter(results, wd != -999) # don't need calms to plot polygons
# add widths
results <- full_join(results, intervals, by = "interval")
results <- mutate(results, n = ifelse(is.na(n), 0, n))
# these results are used for plotting
results2 <- group_by_(results, .dots = c(type, "wd")) %>%
do(calc_poly_coord(., off.set))
angles <- seq(0, 2 * pi, length = 360)
# data for underlying background
datapoly <- data.frame(x = (max.freq + off.set) * sin(angles),
y = (max.freq + off.set) * cos(angles))
# data for gridlines
lines_at <- seq(off.set, mymax, by = myby)
datalines <- data.frame(value = rep(lines_at, each = 360),
x = rep(lines_at, each = 360) * sin(angles),
y = rep(lines_at, each = 360) * cos(angles))
# data for axes
upper <- max.freq + off.set
data.axes <- data.frame(x1 = c(-upper, 0), y1 = c(0, -upper),
x2 = c(upper, 0), y2 = c(0, upper))
plt <- ggplot(results2, aes(x, y, group = interaction(interval, wd))) +
theme_void() +
coord_cartesian(xlim = 1.03 * c(-max.freq - off.set, max.freq + off.set),
ylim = 1.03 * c(-max.freq - off.set, max.freq + off.set)) +
geom_polygon(data = datapoly, aes(x = x, y = y), color = "grey95",
fill = "grey95", inherit.aes = FALSE) +
geom_path(data = datalines, aes(x = x, y = y), color = "white",
inherit.aes = FALSE) +
geom_segment(data = data.axes, aes(x = x1, y = y1, xend = x2, yend = y2),
color = "white", inherit.aes = FALSE) +
annotate("text", upper * -1 * 0.95, 0, label = "W") +
annotate("text", 0, upper * -1 * 0.95, label = "S") +
annotate("text", 0, upper * 0.95, label = "N") +
annotate("text", upper * 0.95, 0, label = "E") +
geom_polygon(aes(fill = interval)) +
scale_fill_manual(values = openColours(cols, max(results$n)),
labels = labs) +
theme(aspect.ratio = 1, strip.background = element_rect(fill = "grey"))
if (!"default" %in% type) {
if (length(type) == 1L)
plt <- plt + facet_wrap(reformulate(type), labeller = label_parsed)
if (length(type) == 2L)
plt <- plt + facet_grid(paste(type[2], "~", type[1]), labeller = label_parsed)
}
plot(plt)
newdata <- results
output <- list(plot = plt, data = newdata, call = match.call())
class(output) <- "openair"
invisible(output)
}
## adds a line showing probability wind direction is from a particular sector
## used when normalise = TRUE
panel.wdprob <- function(dat, seg, angle, off.set) {
len1 <- off.set
x.off <- 0; y.off <- 0
makeline <- function(i, dat) {
theta <- seq((dat$wd[i] - seg * angle / 2), (dat$wd[i] + seg * angle / 2),
length.out = (angle - 2) * 10)
theta <- ifelse(theta < 1, 360 - theta, theta)
theta <- theta * pi / 180
x1 <- len1 * sin(theta) + x.off
x2 <- rev((dat$norm[i] + off.set) * sin(theta) + x.off)
y1 <- len1 * cos(theta) + x.off
y2 <- rev((dat$norm[i] + off.set) * cos(theta) + x.off)
lpolygon(c(x1, x2), c(y1, y2), col = "transparent", border = "black", lwd = 2)
}
lapply(1:nrow(dat), makeline, dat)
}
calc_poly_coord <- function(data, off.set) {
# order by n
data <- arrange(data, n)
data$wd <- data$wd * pi / 180
data$value <- data$value + off.set
# groups of 4 points to make a polygon, then back to starting point
coord_x <- sapply(2:(nrow(data)),
function (i)
with(data,
c(value[i - 1] * sin(wd[i - 1]) -
width[i] * cos(wd[i - 1]),
value[i - 1] * sin(wd[i - 1]) +
width[i] * cos(wd[i - 1]),
value[i] * sin(wd[i]) +
width[i] * cos(wd[i]),
value[i] * sin(wd[i]) -
width[i] * cos(wd[i]),
value[i - 1] * sin(wd[i - 1]) -
width[i] * cos(wd[i - 1]))))
coord_x <- as.vector(coord_x)
coord_y <- sapply(2:(nrow(data)),
function (i)
with(data,
c(value[i - 1] * cos(wd[i - 1]) +
width[i] * sin(wd[i - 1]),
value[i - 1] * cos(wd[i - 1]) -
width[i] * sin(wd[i - 1]),
value[i] * cos(wd[i]) -
width[i] * sin(wd[i]),
value[i] * cos(wd[i]) +
width[i] * sin(wd[i]),
value[i - 1] * cos(wd[i - 1]) +
width[i] * sin(wd[i - 1]))))
coord_y <- as.vector(coord_y)
# make a data frame with all the coordinates for the polygons
result <- data_frame(x = coord_x,
y = coord_y,
interval = rep(data$interval[2:nrow(data)],
each = 5))
return(result)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.