R/wind_rose_plot_routine.R

In MariekeDirk/WINS50_Lidar: Lidar data analysis for the WINS50 project

#' Wind rose plotting routine
#' @description Based on ggplot2 plotting routine
#' @param spd numerical values of the speed
#' @param dir direction
#' @param spdres slipt number
#' @param dirres resolution of the number of directions
#' @param spdmin minimum speed
#' @param spdmax maximum speed
#' @param spdseq NULL
#' @param palette color palette to be used
#' @param countmax NA
#' @param debug to debug the code
#' @importFrom rlang .data
#' @author Marieke Dirksen
#' @export
plot_windrose <- function(spd,
                          dir,
                          spdres = 2,
                          dirres = 30,
                          spdmin = 2,
                          spdmax = 20,
                          spdseq = NULL,
                          palette = "YlGnBu",
                          countmax = NA,
                          debug = 0){


  if(length(spd)!=length(dir)){
    message("spd and dir have different lengths, returning FALSE")
    return(FALSE)
  }
  # 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{
    message("spd or dir are non numeric input vectors, returning FALSE")
    return(FALSE)
  }

  # 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 <- grDevices::colorRampPalette(RColorBrewer::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. <- stats::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(utils::packageVersion("ggplot2") > "2.2"){
    cat("Hadley broke my code\n")
    data$spd.binned = with(data, factor(spd.binned, levels = rev(levels(spd.binned))))#.data$
    spd.colors = rev(spd.colors)
  }

  # create the plot ----
  p.windrose <- ggplot2::ggplot(data = data,
                       ggplot2::aes(x = dir.binned,
                           fill = spd.binned)) + #.data$
    ggplot2::geom_bar() +
    ggplot2::theme_bw()+
    ggplot2::scale_x_discrete(drop = FALSE,
                     labels = ggplot2::waiver()) +
    ggplot2::coord_polar(start = -((dirres/2)/360) * 2*pi) +
    ggplot2::scale_fill_manual(name = "Wind Speed (m/s)",
                      values = spd.colors,
                      drop = FALSE) +
    ggplot2::theme(axis.title.x = ggplot2::element_blank(),
          #panel.border = element_rect(colour = "blank"),
          panel.grid.major = ggplot2::element_line(colour="grey65"))

  # adjust axes if required
  if (!is.na(countmax)){
    p.windrose <- p.windrose +
      ggplot2::ylim(c(0,countmax))
  }

  # print the plot
  print(p.windrose)

  # return the handle to the wind rose
  return(p.windrose)
}