为了论证，我们假设我们正在使用data.in 数据框，它有两个数据列和某种日期/时间信息。我们最初会忽略日期和时间信息。
ggplot 函数
我已经编写了下面的函数。我对其他人的经验或关于如何改进这一点的建议感兴趣。
# WindRose.R
require(ggplot2)
require(RColorBrewer)
plot.windrose <- function(data,
spd,
dir,
spdres = 2,
dirres = 30,
spdmin = 2,
spdmax = 20,
spdseq = NULL,
palette = "YlGnBu",
countmax = NA,
debug = 0){
# Look to see what data was passed in to the function
if (is.numeric(spd) & is.numeric(dir)){
# assume that we've been given vectors of the speed and direction vectors
data <- data.frame(spd = spd,
dir = dir)
spd = "spd"
dir = "dir"
} else if (exists("data")){
# Assume that we've been given a data frame, and the name of the speed
# and direction columns. This is the format we want for later use.
}
# Tidy up input data ----
n.in <- NROW(data)
dnu <- (is.na(data[[spd]]) | is.na(data[[dir]]))
data[[spd]][dnu] <- NA
data[[dir]][dnu] <- NA
# figure out the wind speed bins ----
if (missing(spdseq)){
spdseq <- seq(spdmin,spdmax,spdres)
} else {
if (debug >0){
cat("Using custom speed bins \n")
}
}
# get some information about the number of bins, etc.
n.spd.seq <- length(spdseq)
n.colors.in.range <- n.spd.seq - 1
# create the color map
spd.colors <- colorRampPalette(brewer.pal(min(max(3,
n.colors.in.range),
min(9,
n.colors.in.range)),
palette))(n.colors.in.range)
if (max(data[[spd]],na.rm = TRUE) > spdmax){
spd.breaks <- c(spdseq,
max(data[[spd]],na.rm = TRUE))
spd.labels <- c(paste(c(spdseq[1:n.spd.seq-1]),
'-',
c(spdseq[2:n.spd.seq])),
paste(spdmax,
"-",
max(data[[spd]],na.rm = TRUE)))
spd.colors <- c(spd.colors, "grey50")
} else{
spd.breaks <- spdseq
spd.labels <- paste(c(spdseq[1:n.spd.seq-1]),
'-',
c(spdseq[2:n.spd.seq]))
}
data$spd.binned <- cut(x = data[[spd]],
breaks = spd.breaks,
labels = spd.labels,
ordered_result = TRUE)
# clean up the data
data. <- na.omit(data)
# figure out the wind direction bins
dir.breaks <- c(-dirres/2,
seq(dirres/2, 360-dirres/2, by = dirres),
360+dirres/2)
dir.labels <- c(paste(360-dirres/2,"-",dirres/2),
paste(seq(dirres/2, 360-3*dirres/2, by = dirres),
"-",
seq(3*dirres/2, 360-dirres/2, by = dirres)),
paste(360-dirres/2,"-",dirres/2))
# assign each wind direction to a bin
dir.binned <- cut(data[[dir]],
breaks = dir.breaks,
ordered_result = TRUE)
levels(dir.binned) <- dir.labels
data$dir.binned <- dir.binned
# Run debug if required ----
if (debug>0){
cat(dir.breaks,"\n")
cat(dir.labels,"\n")
cat(levels(dir.binned),"\n")
}
# deal with change in ordering introduced somewhere around version 2.2
if(packageVersion("ggplot2") > "2.2"){
cat("Hadley broke my code\n")
data$spd.binned = with(data, factor(spd.binned, levels = rev(levels(spd.binned))))
spd.colors = rev(spd.colors)
}
# create the plot ----
p.windrose <- ggplot(data = data,
aes(x = dir.binned,
fill = spd.binned)) +
geom_bar() +
scale_x_discrete(drop = FALSE,
labels = waiver()) +
coord_polar(start = -((dirres/2)/360) * 2*pi) +
scale_fill_manual(name = "Wind Speed (m/s)",
values = spd.colors,
drop = FALSE) +
theme(axis.title.x = element_blank())
# adjust axes if required
if (!is.na(countmax)){
p.windrose <- p.windrose +
ylim(c(0,countmax))
}
# print the plot
print(p.windrose)
# return the handle to the wind rose
return(p.windrose)
}
概念和逻辑证明
我们现在将检查代码是否符合我们的预期。为此，我们将使用简单的演示数据集。
# try the default settings
p0 <- plot.windrose(spd = data.in$ws.80,
dir = data.in$wd.80)
这给了我们这个情节：
所以：我们已经按照风向和风速正确地对数据进行了分类，并且按照预期对超出范围的数据进行了编码。看起来不错！
使用此功能
现在我们加载真实数据。我们可以从 URL 加载它：
data.in <- read.csv(file = "http://midcdmz.nrel.gov/apps/plot.pl?site=NWTC&start=20010824&edy=26&emo=3&eyr=2062&year=2013&month=1&day=1&endyear=2013&endmonth=12&endday=31&time=0&inst=21&inst=39&type=data&wrlevel=2&preset=0&first=3&math=0&second=-1&value=0.0&user=0&axis=1",
col.names = c("date","hr","ws.80","wd.80"))
或来自文件：
data.in <- read.csv(file = "A:/blah/20130101.csv",
col.names = c("date","hr","ws.80","wd.80"))
快捷方式
将它与 M2 数据一起使用的简单方法是为 spd 和 dir（速度和方向）传入单独的向量：
# try the default settings
p1 <- plot.windrose(spd = data.in$ws.80,
dir = data.in$wd.80)
这给了我们这个情节：
如果我们想要自定义垃圾箱，我们可以将它们添加为参数：
p2 <- plot.windrose(spd = data.in$ws.80,
dir = data.in$wd.80,
spdseq = c(0,3,6,12,20))
使用数据框和列名
为了使绘图更兼容ggplot()，您还可以传入一个数据框以及速度和方向变量的名称：
p.wr2 <- plot.windrose(data = data.in,
spd = "ws.80",
dir = "wd.80")
由另一个变量分面
我们还可以使用 ggplot 的分面功能按月或按年绘制数据。让我们从data.in 中的日期和小时信息中获取时间戳开始，然后转换为月份和年份：
# first create a true POSIXCT timestamp from the date and hour columns
data.in$timestamp <- as.POSIXct(paste0(data.in$date, " ", data.in$hr),
tz = "GMT",
format = "%m/%d/%Y %H:%M")
# Convert the time stamp to years and months
data.in$Year <- as.numeric(format(data.in$timestamp, "%Y"))
data.in$month <- factor(format(data.in$timestamp, "%B"),
levels = month.name)
然后你可以应用分面来显示风玫瑰如何随月变化：
# recreate p.wr2, so that includes the new data
p.wr2 <- plot.windrose(data = data.in,
spd = "ws.80",
dir = "wd.80")
# now generate the faceting
p.wr3 <- p.wr2 + facet_wrap(~month,
ncol = 3)
# and remove labels for clarity
p.wr3 <- p.wr3 + theme(axis.text.x = element_blank(),
axis.title.x = element_blank())
评论
关于该功能及其使用方法的一些注意事项：
- 输入是：
- 速度向量 (
spd) 和方向 (dir) 或数据框的名称以及包含速度和方向数据的列的名称。
- 风速 (
spdres) 和风向 (dirres) 的 bin 大小可选值。
-
palette 是 colorbrewer 顺序调色板的名称，
-
countmax 设置风玫瑰的范围。
-
debug 是一个开关 (0,1,2)，用于启用不同级别的调试。
- 我希望能够设置绘图的最大速度 (
spdmax) 和计数 (countmax)，以便比较来自不同数据集的风玫瑰
- 如果风速超过 (
spdmax)，则会将其添加为灰色区域（见图）。我可能也应该编写类似spdmin 的代码，并对风速小于该值的区域进行颜色代码。
- 根据请求，我实施了一种使用自定义风速箱的方法。可以使用
spdseq = c(1,3,5,12) 参数添加它们。
- 您可以使用常用的 ggplot 命令删除度数 bin 标签以清除 x 轴：
p.wr3 + theme(axis.text.x = element_blank(),axis.title.x = element_blank())。
- 最近 ggplot2 更改了 bin 的顺序，因此绘图无法正常工作。我认为这是 2.2 版。但是，如果您的绘图看起来有点奇怪，请更改代码，以便“2.2”的测试可能是“2.1”或“2.0”。