diff --git a/NEWS.md b/NEWS.md
index 92f55549..7b036f76 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -24,6 +24,8 @@
   
   - `source = "all"` will return all available metadata (including KCL and Europe).
 
+- new function `runRegression()` for extracting 'dilution lines' from air quality and other data. Online manual will be updated with principles and examples. 
+
 - `calendarPlot()` now automatically creates its own `labels` if `breaks` are specified. For example, `c(0, 10, 20)` will create the labels `c("0 - 10", "10 - 20")`. `labels` can still be used to override the default values. (#341)
 
 - return tibble from `timeAverage()`.
@@ -32,9 +34,7 @@
 
 - Move regression formula off main plot for `polarPlot()` for clarity and label slope as 'm'.
 
-- new function `runRegression()` for extracting 'dilution lines' from air quality and other data. Online manual will be updated with principles and examples. 
-
-- Tweak seasonal trend decomposition using STL to allow the seasonal amplitude to vary more. Affects `smoothTrend` and `TheilSen`.
+- Tweak seasonal trend decomposition using STL to allow the seasonal amplitude to vary more. Affects `smoothTrend()` and `TheilSen()`.
 
 ## Bug Fixes
 
@@ -44,6 +44,8 @@
 
 - The `year` argument of `importMeta()` is now respected when `source = "kcl"` and `"europe"`.
 
+- Several of the directional analysis plot family (e.g., `polarFreq()`) have been refactored to use `is.null()` or `is.na()` over `missing()`. While predominantly an internal change, this should be make these functions easier to use inside of other functions (e.g., `function(data, breaks = NA) polarFreq(data, breaks = breaks))` will now run successfully).
+
 - For `calendarPlot` when annotated with ws or wd arrows, use max ws/wd that corresponds to hour of maximum pollutant concentration and not simple the max ws/wd for a day.
 
 # openair 2.17-0
diff --git a/R/percentileRose.R b/R/percentileRose.R
index 1f361830..8f990adc 100644
--- a/R/percentileRose.R
+++ b/R/percentileRose.R
@@ -453,9 +453,7 @@ percentileRose <- function(
 
   ## nice intervals for pollutant concentrations
   tmp <- (results.grid$x ^ 2 + results.grid$y ^ 2) ^ 0.5
-  if (missing(intervals)) intervals <- pretty(c(min(tmp, na.rm = TRUE), max(tmp, na.rm = TRUE)))
-
-
+  if (is.null(intervals)) intervals <- pretty(c(min(tmp, na.rm = TRUE), max(tmp, na.rm = TRUE)))
 
   labs <- intervals ## the labels
 
diff --git a/R/polarAnnulus.R b/R/polarAnnulus.R
index af1d9cd1..99041d77 100644
--- a/R/polarAnnulus.R
+++ b/R/polarAnnulus.R
@@ -150,7 +150,7 @@ polarAnnulus <-
            type = "default",
            statistic = "mean",
            percentile = NA,
-           limits = c(0, 100),
+           limits = NULL,
            cols = "default",
            width = "normal",
            min.bin = 1,
@@ -159,7 +159,7 @@ polarAnnulus <-
            force.positive = TRUE,
            k = c(20, 10),
            normalise = FALSE,
-           key.header = "",
+           key.header = statistic,
            key.footer = pollutant,
            key.position = "right",
            key = TRUE,
@@ -183,7 +183,6 @@ polarAnnulus <-
     percentile <- 50
   }
 
-  if (missing(key.header)) key.header <- statistic
   if (key.header[1] == "weighted.mean") key.header <- c("weighted", "mean")
   if (key.header[1] == "percentile") key.header <- c(paste(percentile, "th", sep = ""), "percentile")
   if (key.header[1] == "cpf") key.header <- c("CPF", "probability")
@@ -510,7 +509,7 @@ polarAnnulus <-
   ## auto-scaling
   nlev <- 200 # preferred number of intervals
   ## handle missing breaks arguments
-  if (missing(limits)) {
+  if (any(is.null(limits)) | any(is.na(limits))) {
     # breaks <- pretty(results.grid$z, n = nlev)
     breaks <- seq(min(results.grid$z, na.rm = TRUE), max(results.grid$z, na.rm = TRUE),
                   length.out = nlev
diff --git a/R/polarDiff.R b/R/polarDiff.R
index 7c1fce49..d9611070 100644
--- a/R/polarDiff.R
+++ b/R/polarDiff.R
@@ -42,10 +42,13 @@ polarDiff <- function(
     after,
     pollutant = "nox",
     x = "ws",
-    limits = NA,
+    limits = NULL,
     plot = TRUE,
     ...
 ) {
+  if (is.null(limits)){
+    limits <- NA
+  }
 
   # extra args setup
   Args <- list(...)
diff --git a/R/polarFreq.R b/R/polarFreq.R
index e1c6088d..2b16f76c 100644
--- a/R/polarFreq.R
+++ b/R/polarFreq.R
@@ -116,12 +116,12 @@
 #' \dontrun{polarFreq(mydata, pollutant = "pm25", statistic
 #' ="weighted.mean", offset = 50, ws.int = 25, trans = FALSE) }
 polarFreq <- function(mydata,
-                      pollutant = "",
+                      pollutant = NULL,
                       statistic = "frequency",
                       ws.int = 1,
                       wd.nint = 36,
                       grid.line = 5,
-                      breaks = seq(0, 5000, 500),
+                      breaks = NULL,
                       cols = "default",
                       trans = TRUE,
                       type = "default",
@@ -183,7 +183,7 @@ polarFreq <- function(mydata,
     trellis.par.set(fontsize = list(text = extra.args$fontsize))
   }
 
-  if (!missing(pollutant)) vars <- c(vars, pollutant)
+  if (!is.null(pollutant)) vars <- c(vars, pollutant)
 
   ## data checks
   mydata <- checkPrep(mydata, vars, type, remove.calm = FALSE)
@@ -198,19 +198,20 @@ polarFreq <- function(mydata,
   mydata <- cutData(mydata, type, ...)
 
   ## if pollutant chosen but no statistic - use mean, issue warning
-  if (!missing(pollutant) & missing(statistic)) {
+  if (statistic == "frequency" & !is.null(pollutant)) {
+    cli::cli_warn(c("x" = "{.code statistic == 'frequency'} incompatible with a defined {.field pollutant}.",
+                    "i" = "Setting {.field statistic} to {.code 'mean'}."))
     statistic <- "mean"
-    warning("No statistic chosen, using mean")
   }
 
   ## if statistic chosen but no pollutant stop
-  if (!missing(statistic) & missing(pollutant)) {
-    stop("No pollutant chosen, please choose one e.g. pollutant = 'nox'")
+  if (statistic != "frequency" & is.null(pollutant)) {
+    cli::cli_abort(c("x" = "No {.field pollutant} chosen",
+                     "i" = "Please choose a {.field pollutant}, e.g., {.code pollutant = 'nox'}"))
   }
 
-  if (!missing(breaks)) trans <- FALSE ## over-ride transform if breaks supplied
+  if (!(any(is.null(breaks)) | any(is.na(breaks)))) trans <- FALSE ## over-ride transform if breaks supplied
 
-  if (missing(key.header)) key.header <- statistic
   if (key.header == "weighted.mean") key.header <- c("contribution", "(%)")
 
   ## apply square root transform?
@@ -319,7 +320,7 @@ polarFreq <- function(mydata,
 
   nlev <- 200
   ## handle missing breaks arguments
-  if (missing(breaks)) {
+  if (any(is.null(breaks)) | any(is.na(breaks))) {
     breaks <- unique(c(0, pretty(results.grid$weights, nlev)))
     br <- pretty((c(0, results.grid$weights) ^ coef), n = 10) ## breaks for scale
   } else {
diff --git a/R/polarPlot.R b/R/polarPlot.R
index da3e25a9..885a2792 100644
--- a/R/polarPlot.R
+++ b/R/polarPlot.R
@@ -323,9 +323,9 @@
 #' @param kernel Type of kernel used for the weighting procedure for when
 #'   correlation or regression techniques are used. Only \code{"gaussian"} is
 #'   supported but this may be enhanced in the future.
-#'   
+#'
 #' @param formula.label When pair-wise statistics such as regression slopes are
-#'   calculated and plotted, should a formula label be displayed? 
+#'   calculated and plotted, should a formula label be displayed?
 #'
 #' @param tau The quantile to be estimated when \code{statistic} is set to
 #'   \code{"quantile.slope"}. Default is \code{0.5} which is equal to the median
@@ -432,7 +432,7 @@ polarPlot <-
            wd = "wd",
            type = "default",
            statistic = "mean",
-           limits = NA,
+           limits = NULL,
            exclude.missing = TRUE,
            uncertainty = FALSE,
            percentile = NA,
@@ -446,7 +446,7 @@ polarPlot <-
            force.positive = TRUE,
            k = 100,
            normalise = FALSE,
-           key.header = "",
+           key.header = statistic,
            key.footer = pollutant,
            key.position = "right",
            key = TRUE,
@@ -461,7 +461,6 @@ polarPlot <-
            alpha = 1,
            plot = TRUE,
            ...) {
-
     ## get rid of R check annoyances
     z <- . <- NULL
 
@@ -490,8 +489,9 @@ polarPlot <-
     }
 
     # if statistic is trend, then don't force to be positive
-    if (statistic == "trend")
+    if (statistic == "trend") {
       force.positive <- FALSE
+    }
 
     # names of variables for use later
     nam.x <- x
@@ -515,7 +515,6 @@ polarPlot <-
 
     if (length(weights) != 3) stop("weights should be of length 3.")
 
-    if (missing(key.header)) key.header <- statistic
     if (key.header == "nwr") key.header <- "NWR"
     if (key.header == "weighted_mean") key.header <- "weighted\nmean"
     if (key.header == "percentile") {
@@ -579,8 +578,9 @@ polarPlot <-
     ## extract variables of interest
     vars <- c(wd, x, pollutant)
 
-    if (statistic == "york_slope")
+    if (statistic == "york_slope") {
       vars <- c(vars, x_error, y_error)
+    }
 
     if (any(type %in% dateTypes)) vars <- c(vars, "date")
 
@@ -639,7 +639,7 @@ polarPlot <-
     min.ws <- min(mydata[[x]], na.rm = TRUE)
     clip <- TRUE ## used for removing data where ws > upper
 
-    if (missing(upper)) {
+    if (is.na(upper)) {
       upper <- max.ws
       clip <- FALSE
     }
@@ -675,7 +675,6 @@ polarPlot <-
     )
 
     if (statistic == "cpf") {
-
       ## can be interval of percentiles or a single (threshold)
       if (length(percentile) > 1) {
         statistic <- "cpfi" # CPF interval
@@ -725,7 +724,6 @@ polarPlot <-
 
 
     prepare.grid <- function(mydata) {
-
       ## identify which ws and wd bins the data belong
       wd <- cut(
         wd.int * ceiling(mydata[[wd]] / wd.int - 0.5),
@@ -739,8 +737,7 @@ polarPlot <-
       )
 
       if (!statistic %in% c(correlation_stats, "nwr", "trend")) {
-        binned <- switch(
-          statistic,
+        binned <- switch(statistic,
           frequency = tapply(mydata[[pollutant]], list(wd, x), function(x) {
             length(na.omit(x))
           }),
@@ -794,7 +791,6 @@ polarPlot <-
           ))
 
         binned <- binned$conc
-
       } else if (toupper(statistic) == "TREND") {
         binned <- rowwise(ws.wd) %>%
           summarise(simple_kernel_trend(
@@ -806,9 +802,7 @@ polarPlot <-
           ))
 
         binned <- binned$conc
-
       } else {
-
         binned <- rowwise(ws.wd) %>%
           summarise(calculate_weighted_statistics(
             across(.cols = everything()),
@@ -851,37 +845,28 @@ polarPlot <-
 
       ## no uncertainty to calculate
       if (!uncertainty) {
-
         ## catch errors when not enough data to calculate surface
         Mgam <- try(gam(binned^n ~ s(u, v, k = k), weights = W), TRUE)
 
         if (!inherits(Mgam, "try-error")) {
-
           pred <- predict.gam(Mgam, input.data)
           pred <- pred^(1 / n)
           pred <- as.vector(pred)
 
           # interpolate results for speed, but not for clustering
           if (extra.args$cluster) {
-
             results <- interp_grid(input.data, z = pred, n = 101)
             int <- 101
-
           } else {
-
             results <- interp_grid(input.data, z = pred, n = 201)
             int <- 201
-
           }
-
-
         } else {
           results <- data.frame(u = u, v = v, z = binned)
           exclude.missing <- FALSE
           warning(call. = FALSE, paste("Not enough data to fit surface.\nTry reducing the value of the smoothing parameter, k to less than ", k, ". \nOr use statistic = 'nwr'.", sep = ""))
         }
       } else {
-
         ## uncertainties calculated, weighted by number of points in each bin
         Mgam <- gam(binned^n ~ s(u, v, k = k), weights = binned.len)
         pred <- predict.gam(Mgam, input.data, se.fit = TRUE)
@@ -898,20 +883,18 @@ polarPlot <-
 
         n <- length(pred)
 
-        # interpolate each uncertainty surface 
-        
-         lower_uncer = interp_grid(input.data, z = Lower, n = 201) %>% mutate(uncertainty = "lower uncertainty")
-         upper_uncer = interp_grid(input.data, z = Upper, n = 201) %>% mutate(uncertainty = "upper uncertainty")
-         prediction = interp_grid(input.data, z = pred, n = 201) %>% mutate(uncertainty = "prediction")
-         results <- bind_rows(prediction, lower_uncer, upper_uncer)
-         int <- 201
-        
+        # interpolate each uncertainty surface
+
+        lower_uncer <- interp_grid(input.data, z = Lower, n = 201) %>% mutate(uncertainty = "lower uncertainty")
+        upper_uncer <- interp_grid(input.data, z = Upper, n = 201) %>% mutate(uncertainty = "upper uncertainty")
+        prediction <- interp_grid(input.data, z = pred, n = 201) %>% mutate(uncertainty = "prediction")
+        results <- bind_rows(prediction, lower_uncer, upper_uncer)
+        int <- 201
       }
 
 
       ## function to remove points too far from original data
       exclude <- function(results) {
-
         ## exclude predictions too far from data (from mgcv)
         x <- seq(-upper, upper, length = int)
         y <- x
@@ -993,10 +976,11 @@ polarPlot <-
       id <- which(res$z < -1)
       if (length(id) > 0) res$z[id] <- -1
     }
-    
+
     # if regression, set key.footer to 'm' (slope)
-    if (grepl("slope|intercept", statistic) & length(pollutant == 2)) 
+    if (grepl("slope|intercept", statistic) & length(pollutant == 2)) {
       key.footer <- "m"
+    }
 
 
     # Labels for correlation and regression, keep lower case like other labels
@@ -1023,8 +1007,7 @@ polarPlot <-
 
     ## handle missing breaks arguments
 
-    if (is.na(limits[1])) {
-
+    if (any(is.null(limits)) | any(is.na(limits))) {
       # breaks <- pretty(res$z, n = nlev)
       breaks <- seq(
         min(res$z, na.rm = TRUE), max(res$z, na.rm = TRUE),
@@ -1034,7 +1017,6 @@ polarPlot <-
       labs <- labs[labs >= min(breaks) & labs <= max(breaks)]
       at <- labs
     } else {
-
       ## handle user limits and clipping
       breaks <- seq(min(limits), max(limits), length.out = nlev)
       labs <- pretty(breaks, 7)
@@ -1092,10 +1074,11 @@ polarPlot <-
       labels <- labels[-1]
       intervals <- intervals[-1]
     }
-    
+
     # if uncertainty = TRUE, change type for plotting (3 panels)
-    if (uncertainty)
+    if (uncertainty) {
       type <- "uncertainty"
+    }
 
     temp <- paste(type, collapse = "+")
     myform <- formula(paste("z ~ u * v | ", temp, sep = ""))
@@ -1115,24 +1098,24 @@ polarPlot <-
       ylim = c(-upper * 1.025, upper * 1.025),
       colorkey = FALSE,
       legend = legend,
-      
+
       # add footnote formula if regression used
-      page = function(n)
-        if (formula.label & grepl("slope|intercept", statistic) & 
-            length(pollutant == 2)) {
-        grid.text(
-          quickText(paste0(
-            "Formula: ", pollutant[1], " = m.", pollutant[2], " + c"
-          )),
-          x = .99,
-          y = 0.01,
-          default.units = "npc",
-          gp = gpar(fontsize = 10),
-          just = c("right", "bottom")
-        )},
-      
+      page = function(n) {
+        if (formula.label & grepl("slope|intercept", statistic) &
+          length(pollutant == 2)) {
+          grid.text(
+            quickText(paste0(
+              "Formula: ", pollutant[1], " = m.", pollutant[2], " + c"
+            )),
+            x = .99,
+            y = 0.01,
+            default.units = "npc",
+            gp = gpar(fontsize = 10),
+            just = c("right", "bottom")
+          )
+        }
+      },
       panel = function(x, y, z, subscripts, ...) {
-
         ## show missing data due to min.bin
         if (min.bin > 1) {
           panel.levelplot(
@@ -1183,7 +1166,6 @@ polarPlot <-
         ltext(0.07 * upper, upper * -1 * 0.95, "S", cex = 0.7)
         ltext(0.07 * upper, upper * 0.95, "N", cex = 0.7)
         ltext(upper * 0.95, 0.07 * upper, "E", cex = 0.7)
-
       }
     )
 
@@ -1198,8 +1180,9 @@ polarPlot <-
         plot(plt)
       } else {
         plot(useOuterStrips(plt,
-                            strip = strip,
-                            strip.left = strip.left))
+          strip = strip,
+          strip.left = strip.left
+        ))
       }
     }
 
@@ -1213,10 +1196,8 @@ polarPlot <-
 
 # Gaussian bivariate density function
 gauss_dens <- function(x, y, mx, my, sx, sy) {
-  
-  (1 / (2 * pi * sx *sy )) *
-    exp((-1/2) * ((x - mx) ^ 2 / sx ^ 2 + (y - my) ^ 2 / sy^2))
-  
+  (1 / (2 * pi * sx * sy)) *
+    exp((-1 / 2) * ((x - mx)^2 / sx^2 + (y - my)^2 / sy^2))
 }
 
 # NWR kernel calculations
@@ -1226,40 +1207,38 @@ simple_kernel <- function(data, mydata, x = "ws",
   # Centres
   ws1 <- data[[1]]
   wd1 <- data[[2]]
-  
+
   # centred ws, wd
   ws_cent <- mydata[[x]] - ws1
   wd_cent <- mydata[[y]] - wd1
-  wd_cent = ifelse(wd_cent < -180, wd_cent + 360, wd_cent)
-  
+  wd_cent <- ifelse(wd_cent < -180, wd_cent + 360, wd_cent)
+
   weight <- gauss_dens(ws_cent, wd_cent, 0, 0, ws_spread, wd_spread)
 
   conc <- sum(mydata[[pollutant]] * weight) /
     sum(weight)
-  
+
   return(data.frame(conc = conc))
 }
 
 # function to to kernel weighting of a quantile regression trend
 simple_kernel_trend <- function(data, mydata, x = "ws",
-                          y = "wd", pollutant, date,
-                          ws_spread, wd_spread, kernel, tau) {
+                                y = "wd", pollutant, date,
+                                ws_spread, wd_spread, kernel, tau) {
   # Centres
   ws1 <- data[[1]]
   wd1 <- data[[2]]
 
   # Gaussian bivariate density function
   gauss_dens <- function(x, y, mx, my, sx, sy) {
-
-    (1 / (2 * pi * sx *sy )) *
-      exp((-1/2) * ((x - mx) ^ 2 / sx ^ 2 + (y - my) ^ 2 / sy^2))
-
+    (1 / (2 * pi * sx * sy)) *
+      exp((-1 / 2) * ((x - mx)^2 / sx^2 + (y - my)^2 / sy^2))
   }
 
   # centred ws, wd
   ws_cent <- mydata[[x]] - ws1
   wd_cent <- mydata[[y]] - wd1
-  wd_cent = ifelse(wd_cent < -180, wd_cent + 360, wd_cent)
+  wd_cent <- ifelse(wd_cent < -180, wd_cent + 360, wd_cent)
 
   weight <- gauss_dens(ws_cent, wd_cent, 0, 0, ws_spread, wd_spread)
   weight <- weight / max(weight)
@@ -1286,12 +1265,10 @@ simple_kernel_trend <- function(data, mydata, x = "ws",
 
   # Extract statistics
   if (!inherits(fit, "try-error")) {
-
     # Extract statistics
-    slope = 365.25 * 24 * 3600 * fit$coefficients[2]
-
+    slope <- 365.25 * 24 * 3600 * fit$coefficients[2]
   } else {
-    slope <-  NA
+    slope <- NA
   }
 
 
@@ -1305,163 +1282,155 @@ calculate_weighted_statistics <-
            y = "wd", pol_1, pol_2,
            ws_spread, wd_spread, kernel, tau,
            x_error, y_error) {
-  weight <- NULL
-  # Centres
-  ws1 <- data[[1]]
-  wd1 <- data[[2]]
-
-  # centred ws, wd
-  ws_cent <- mydata[[x]] - ws1
-  wd_cent <- mydata[[y]] - wd1
-  wd_cent = ifelse(wd_cent < -180, wd_cent + 360, wd_cent)
-
-  weight <- gauss_dens(ws_cent, wd_cent, 0, 0, ws_spread, wd_spread)
-  weight <- weight / max(weight)
-
-  mydata$weight <- weight
-
-  # Select and filter
-  vars <- c(pol_1, pol_2, "weight")
-
-  if (!all(is.na(c(x_error, y_error))))
-    vars <- c(vars, x_error, y_error)
-  thedata <- mydata[vars]
-#  thedata <- thedata[complete.cases(thedata), ]
-
-  # don't fit all data - takes too long with no gain
-  thedata <- thedata[which(thedata$weight > 0.001), ]
+    weight <- NULL
+    # Centres
+    ws1 <- data[[1]]
+    wd1 <- data[[2]]
 
-  # don't try and calculate stats is there's not much data
-  if (nrow(thedata) < 100) return(data.frame(ws1, wd1, NA))
+    # centred ws, wd
+    ws_cent <- mydata[[x]] - ws1
+    wd_cent <- mydata[[y]] - wd1
+    wd_cent <- ifelse(wd_cent < -180, wd_cent + 360, wd_cent)
 
-  # useful for showing what the weighting looks like as a surface
-  # openair::scatterPlot(mydata, x = "ws", y = "wd", z = "weight", method = "level")
+    weight <- gauss_dens(ws_cent, wd_cent, 0, 0, ws_spread, wd_spread)
+    weight <- weight / max(weight)
 
+    mydata$weight <- weight
 
-  if (statistic %in% c("r", "Pearson","Spearman")) {
+    # Select and filter
+    vars <- c(pol_1, pol_2, "weight")
 
-    if (statistic == "r") statistic <- "Pearson"
-
-    # Weighted Pearson correlation
-    stat_weighted <- contCorr(thedata[[pol_1]], thedata[[pol_2]],
-      w = thedata$weight, method = statistic
-    )
-
-    result <- data.frame(ws1, wd1, stat_weighted)
-  }
+    if (!all(is.na(c(x_error, y_error)))) {
+      vars <- c(vars, x_error, y_error)
+    }
+    thedata <- mydata[vars]
+    #  thedata <- thedata[complete.cases(thedata), ]
 
-  # Simple least squared regression with weights
-  if (statistic %in% c("slope", "intercept")) {
+    # don't fit all data - takes too long with no gain
+    thedata <- thedata[which(thedata$weight > 0.001), ]
 
-    # Drop dplyr's data frame for formula
-    thedata <- data.frame(thedata)
+    # don't try and calculate stats is there's not much data
+    if (nrow(thedata) < 100) {
+      return(data.frame(ws1, wd1, NA))
+    }
 
-    # Calculate model, no warnings on perfect fits.
-    fit <- lm(thedata[, pol_1] ~ thedata[, pol_2], weights = thedata[, "weight"])
+    # useful for showing what the weighting looks like as a surface
+    # openair::scatterPlot(mydata, x = "ws", y = "wd", z = "weight", method = "level")
 
-    # Extract statistics
-    if (statistic == "slope") stat_weighted <- fit$coefficients[2]
-    if (statistic == "intercept") stat_weighted <- fit$coefficients[1]
 
-    # Bind together
-    result <- data.frame(ws1, wd1, stat_weighted)
-  }
+    if (statistic %in% c("r", "Pearson", "Spearman")) {
+      if (statistic == "r") statistic <- "Pearson"
 
-  if (statistic == "york_slope") {
+      # Weighted Pearson correlation
+      stat_weighted <- contCorr(thedata[[pol_1]], thedata[[pol_2]],
+        w = thedata$weight, method = statistic
+      )
 
-    thedata <- as.data.frame(thedata) # so formula works
+      result <- data.frame(ws1, wd1, stat_weighted)
+    }
 
-    result <- try(YorkFit(thedata,
-                          X = names(thedata)[2],
-                          Y = names(thedata)[1],
-                          Xstd = x_error, Ystd = y_error,
-                          weight = thedata$weight), TRUE)
+    # Simple least squared regression with weights
+    if (statistic %in% c("slope", "intercept")) {
+      # Drop dplyr's data frame for formula
+      thedata <- data.frame(thedata)
 
-    # Extract statistics
-    if (!inherits(result, "try-error")) {
+      # Calculate model, no warnings on perfect fits.
+      fit <- lm(thedata[, pol_1] ~ thedata[, pol_2], weights = thedata[, "weight"])
 
       # Extract statistics
-      stat_weighted <- result$Slope
-
-    } else {
-
-      stat_weighted <- NA
+      if (statistic == "slope") stat_weighted <- fit$coefficients[2]
+      if (statistic == "intercept") stat_weighted <- fit$coefficients[1]
 
+      # Bind together
+      result <- data.frame(ws1, wd1, stat_weighted)
     }
 
+    if (statistic == "york_slope") {
+      thedata <- as.data.frame(thedata) # so formula works
 
-    result <- data.frame(ws1, wd1, stat_weighted)
-
-  }
+      result <- try(YorkFit(thedata,
+        X = names(thedata)[2],
+        Y = names(thedata)[1],
+        Xstd = x_error, Ystd = y_error,
+        weight = thedata$weight
+      ), TRUE)
 
-  # Robust linear regression with weights
-  if (grepl("robust", statistic, ignore.case = TRUE)) {
+      # Extract statistics
+      if (!inherits(result, "try-error")) {
+        # Extract statistics
+        stat_weighted <- result$Slope
+      } else {
+        stat_weighted <- NA
+      }
 
-    # Drop dplyr's data frame for formula
-    thedata <- data.frame(thedata)
 
-    # Build model, optimal method (MM) cannot use weights
-    fit <- suppressWarnings(
-      try(MASS::rlm(
-        thedata[, pol_1] ~ thedata[, pol_2],
-        weights = thedata[, "weight"], method = "M"
-      ), TRUE)
-    )
+      result <- data.frame(ws1, wd1, stat_weighted)
+    }
 
+    # Robust linear regression with weights
+    if (grepl("robust", statistic, ignore.case = TRUE)) {
+      # Drop dplyr's data frame for formula
+      thedata <- data.frame(thedata)
+
+      # Build model, optimal method (MM) cannot use weights
+      fit <- suppressWarnings(
+        try(MASS::rlm(
+          thedata[, pol_1] ~ thedata[, pol_2],
+          weights = thedata[, "weight"], method = "M"
+        ), TRUE)
+      )
 
-    # Extract statistics
-    if (!inherits(fit, "try-error")) {
 
       # Extract statistics
-      if (statistic == "robust_slope") stat_weighted <- fit$coefficients[2]
-      if (statistic == "robust_intercept") stat_weighted <- fit$coefficients[1]
-    } else {
-      if (statistic == "robust_slope") stat_weighted <- NA
-      if (statistic == "robust_intercept") stat_weighted <- NA
-    }
-
-    # Bind together
-    result <- data.frame(ws1, wd1, stat_weighted)
-  }
+      if (!inherits(fit, "try-error")) {
+        # Extract statistics
+        if (statistic == "robust_slope") stat_weighted <- fit$coefficients[2]
+        if (statistic == "robust_intercept") stat_weighted <- fit$coefficients[1]
+      } else {
+        if (statistic == "robust_slope") stat_weighted <- NA
+        if (statistic == "robust_intercept") stat_weighted <- NA
+      }
 
+      # Bind together
+      result <- data.frame(ws1, wd1, stat_weighted)
+    }
 
 
-  # Quantile regression with weights
-  if (grepl("quantile", statistic, ignore.case = TRUE)) {
 
-    # quantreg is a Suggests package, so make sure it is there
-    try_require("quantreg", "polarPlot")
+    # Quantile regression with weights
+    if (grepl("quantile", statistic, ignore.case = TRUE)) {
+      # quantreg is a Suggests package, so make sure it is there
+      try_require("quantreg", "polarPlot")
 
-    # Drop dplyr's data frame for formula
-    thedata <- data.frame(thedata)
+      # Drop dplyr's data frame for formula
+      thedata <- data.frame(thedata)
 
-    # Build model
-    suppressWarnings(
-      fit <- try(quantreg::rq(
-        thedata[[pol_1]] ~ thedata[[pol_2]],
-        tau = tau,
-        weights = thedata[["weight"]], method = "br"
-      ), TRUE)
-    )
+      # Build model
+      suppressWarnings(
+        fit <- try(quantreg::rq(
+          thedata[[pol_1]] ~ thedata[[pol_2]],
+          tau = tau,
+          weights = thedata[["weight"]], method = "br"
+        ), TRUE)
+      )
 
-    if (!inherits(fit, "try-error")) {
+      if (!inherits(fit, "try-error")) {
+        # Extract statistics
+        if (statistic == "quantile_slope") stat_weighted <- fit$coefficients[2]
+        if (statistic == "quantile_intercept") stat_weighted <- fit$coefficients[1]
+      } else {
+        if (statistic == "quantile_slope") stat_weighted <- NA
+        if (statistic == "quantile_intercept") stat_weighted <- NA
+      }
 
-      # Extract statistics
-      if (statistic == "quantile_slope") stat_weighted <- fit$coefficients[2]
-      if (statistic == "quantile_intercept") stat_weighted <- fit$coefficients[1]
-    } else {
-      if (statistic == "quantile_slope") stat_weighted <- NA
-      if (statistic == "quantile_intercept") stat_weighted <- NA
+      # Bind together
+      result <- data.frame(ws1, wd1, stat_weighted)
     }
 
-    # Bind together
-    result <- data.frame(ws1, wd1, stat_weighted)
+    # Return
+    result
   }
 
-  # Return
-  result
-}
-
 
 # Taken directly from the boot package to save importing
 corr <- function(d, w = rep(1, nrow(d)) / nrow(d)) {
@@ -1596,14 +1565,14 @@ YorkFit <- function(input_data, X = "X", Y = "Y",
   Ystd <- input_data[[Ystd]]
 
   # don't try regression if < 3 points
-  if (sum(!is.na(X)) < 3 || sum(!is.na(Y)) < 3) return()
+  if (sum(!is.na(X)) < 3 || sum(!is.na(Y)) < 3) {
+    return()
+  }
 
   # used in polar plots - Gaussian kernel weighting
   if (!all(is.na(weight))) {
-
     Xstd <- Xstd / weight
     Ystd <- Ystd / weight
-
   }
 
   Xw <- 1 / (Xstd^2) # X weights
@@ -1618,7 +1587,6 @@ YorkFit <- function(input_data, X = "X", Y = "Y",
   n <- 0 # counter for debugging
 
   while (b.diff > tol && n < 100) {
-
     n <- n + 1 # counter to keep a check on convergence
 
     b.old <- b
@@ -1642,10 +1610,13 @@ YorkFit <- function(input_data, X = "X", Y = "Y",
     b <- sumTOP / sumBOT
 
     # zero or problematic data
-    if (anyNA(b, b.old))
-      return(tibble(Intercept = NA, Slope = NA,
-                    Int_error = NA, Slope_error = NA,
-                    OLS_slope = NA))
+    if (anyNA(b, b.old)) {
+      return(tibble(
+        Intercept = NA, Slope = NA,
+        Int_error = NA, Slope_error = NA,
+        OLS_slope = NA
+      ))
+    }
 
     b.diff <- abs(b - b.old)
   }
@@ -1672,9 +1643,11 @@ YorkFit <- function(input_data, X = "X", Y = "Y",
   sumSint <- sum(wSint, na.rm = TRUE)
   wYorkGOF <- c(sumSint / (lgth - 2), sqrt(2 / (lgth - 2))) # GOF (should equal 1 if assumptions are valid), #standard error in GOF
 
-  ans <- tibble(Intercept = a, Slope = b,
-                Int_error = a.err, Slope_error = b.err,
-                OLS_slope = b0)
+  ans <- tibble(
+    Intercept = a, Slope = b,
+    Int_error = a.err, Slope_error = b.err,
+    OLS_slope = b0
+  )
 
   return(ans)
 }
@@ -1682,7 +1655,6 @@ YorkFit <- function(input_data, X = "X", Y = "Y",
 # bi-linear interpolation on a regular grid
 # allows surface predictions at a low resolution to be interpolated, rather than using a GAM
 interp.surface <- function(obj, loc) {
-
   # obj is a surface or image  object like the list for contour, persp or image.
   # loc a matrix of 2 d locations -- new points to evaluate the surface.
   x <- obj$x
@@ -1709,27 +1681,32 @@ interp.surface <- function(obj, loc) {
   # the four corners of the grid box containing the new
   # points.
   return(z[cbind(lx1, ly1)] * (1 - ex) * (1 - ey) + z[cbind(lx1 +
-                                                              1, ly1)] * ex * (1 - ey) + z[cbind(lx1, ly1 + 1)] * (1 -
-                                                                                                                     ex) * ey + z[cbind(lx1 + 1, ly1 + 1)] * ex * ey)
+    1, ly1)] * ex * (1 - ey) + z[cbind(lx1, ly1 + 1)] * (1 -
+    ex) * ey + z[cbind(lx1 + 1, ly1 + 1)] * ex * ey)
 }
 
 # function to do bilinear interpolation given input grid and number of points required
 interp_grid <- function(input.data, x = "u", y = "v", z, n = 201) {
-
   # current number of points
   int <- length(unique(input.data[[x]]))
 
-  obj <- list(x = seq(min(input.data[[x]]), max(input.data[[x]]), length.out = int),
-              y = seq(min(input.data[[y]]), max(input.data[[y]]), length.out = int),
-              z = matrix(z, nrow = int))
+  obj <- list(
+    x = seq(min(input.data[[x]]), max(input.data[[x]]), length.out = int),
+    y = seq(min(input.data[[y]]), max(input.data[[y]]), length.out = int),
+    z = matrix(z, nrow = int)
+  )
 
-  loc <- expand.grid(x = seq(min(input.data[[x]]), max(input.data[[x]]), length.out = n),
-                     y = seq(min(input.data[[y]]), max(input.data[[y]]), length.out = n))
+  loc <- expand.grid(
+    x = seq(min(input.data[[x]]), max(input.data[[x]]), length.out = n),
+    y = seq(min(input.data[[y]]), max(input.data[[y]]), length.out = n)
+  )
 
   res.interp <- interp.surface(obj, loc)
 
-  out <- expand.grid(u = seq(min(input.data[[x]]), max(input.data[[x]]), length.out = n),
-                     v = seq(min(input.data[[y]]), max(input.data[[y]]), length.out = n))
+  out <- expand.grid(
+    u = seq(min(input.data[[x]]), max(input.data[[x]]), length.out = n),
+    v = seq(min(input.data[[y]]), max(input.data[[y]]), length.out = n)
+  )
 
   results <- data.frame(out, z = res.interp)
 }
diff --git a/man/polarAnnulus.Rd b/man/polarAnnulus.Rd
index 675e259b..c5a3aa90 100644
--- a/man/polarAnnulus.Rd
+++ b/man/polarAnnulus.Rd
@@ -13,7 +13,7 @@ polarAnnulus(
   type = "default",
   statistic = "mean",
   percentile = NA,
-  limits = c(0, 100),
+  limits = NULL,
   cols = "default",
   width = "normal",
   min.bin = 1,
@@ -22,7 +22,7 @@ polarAnnulus(
   force.positive = TRUE,
   k = c(20, 10),
   normalise = FALSE,
-  key.header = "",
+  key.header = statistic,
   key.footer = pollutant,
   key.position = "right",
   key = TRUE,
diff --git a/man/polarDiff.Rd b/man/polarDiff.Rd
index dbce85a7..30d71258 100644
--- a/man/polarDiff.Rd
+++ b/man/polarDiff.Rd
@@ -9,7 +9,7 @@ polarDiff(
   after,
   pollutant = "nox",
   x = "ws",
-  limits = NA,
+  limits = NULL,
   plot = TRUE,
   ...
 )
diff --git a/man/polarFreq.Rd b/man/polarFreq.Rd
index 01f1e73c..82eba5ff 100644
--- a/man/polarFreq.Rd
+++ b/man/polarFreq.Rd
@@ -6,12 +6,12 @@
 \usage{
 polarFreq(
   mydata,
-  pollutant = "",
+  pollutant = NULL,
   statistic = "frequency",
   ws.int = 1,
   wd.nint = 36,
   grid.line = 5,
-  breaks = seq(0, 5000, 500),
+  breaks = NULL,
   cols = "default",
   trans = TRUE,
   type = "default",
diff --git a/man/polarPlot.Rd b/man/polarPlot.Rd
index fc5b6abd..d670656d 100644
--- a/man/polarPlot.Rd
+++ b/man/polarPlot.Rd
@@ -11,7 +11,7 @@ polarPlot(
   wd = "wd",
   type = "default",
   statistic = "mean",
-  limits = NA,
+  limits = NULL,
   exclude.missing = TRUE,
   uncertainty = FALSE,
   percentile = NA,
@@ -25,7 +25,7 @@ polarPlot(
   force.positive = TRUE,
   k = 100,
   normalise = FALSE,
-  key.header = "",
+  key.header = statistic,
   key.footer = pollutant,
   key.position = "right",
   key = TRUE,