diff --git a/R/plottingfns.R b/R/plottingfns.R
index eaceaec..3a7a009 100644
--- a/R/plottingfns.R
+++ b/R/plottingfns.R
@@ -51,7 +51,7 @@ draw_world <- function(g = ggplot() + theme_bw() + xlab("") + ylab(""),inc_borde
     }
 
     ## Now return a gg object with the map overlayed
-    g + geom_path(data = worldmap, aes(x=long, y=lat, group=group), colour="black",size=0.1)
+    g + geom_path(data = worldmap, aes(x=long, y=lat, group=group), colour="black",linewidth=0.1)
 }
 
 #' @rdname show_basis
@@ -234,107 +234,107 @@ EmptyTheme <- function() {
 # ---- FUNCTIONS SPECIFIC TO PLOT() -----
 
 #' @rdname plot
-#' @export 
+#' @export
 setMethod("plot", signature(x = "SRE", y = "list"), function(x, y, ...) {
-  
+
   if(!("newdata" %in% names(y))) stop("y should contain an element named 'newdata'")
-  
+
   plot(x, y, ...)
 })
 
 #' @rdname plot
-#' @export 
+#' @export
 setMethod("plot", signature(x = "SRE", y = "STFDF"), function(x, y, ...) .plot_common(x, y, ...))
 
 #' @rdname plot
-#' @export 
+#' @export
 setMethod("plot", signature(x = "SRE", y = "SpatialPointsDataFrame"), function(x, y, ...) .plot_common(x, y, ...))
 
 #' @rdname plot
-#' @export 
+#' @export
 setMethod("plot", signature(x = "SRE", y = "SpatialPixelsDataFrame"), function(x, y, ...) .plot_common(x, y, ...))
 
 #' @rdname plot
-#' @export 
+#' @export
 setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function(x, y, ...) .plot_common(x, y, ...))
 
 
 .plot_common <- function(object, newdata, ...) {
-  
+
   ## Determine columns we need to plot, and which palette to use
   if (object@method == "EM") {
-    
-    ## Here we will just rename the columns produced by the EM side to 
+
+    ## Here we will just rename the columns produced by the EM side to
     ## align with those of the TMB side.
-    newdata$p_mu <- newdata$mu 
+    newdata$p_mu <- newdata$mu
     newdata$RMSPE_mu <- newdata$sd
     column_names <- c("p_mu", "RMSPE_mu")
     palette <- c("Spectral", "BrBG")
-    
+
   } else if (object@method == "TMB") {
-    
+
     ## Unfortunately, we don't know the value of type specified in predict().
-    ## We will just check the existence of quantities using if statements. 
+    ## We will just check the existence of quantities using if statements.
     ## Usually, people want prediction interval widths, but the output of predict
     ## doesn't give interval widths directly. Instead, it gives user-specified
     ## percentiles, so here we construct some interval widths assuming the user
     ## has specified the 5th and 95th percentiles (default).
-    
-    ## If the 5th and 95th percentiles are present, construct the 90% predictive interval width. 
-    if (all(c("Y_percentile_5","Y_percentile_95") %in% names(newdata@data))) 
-      newdata@data$interval90_Y <- newdata$Y_percentile_95 - newdata$Y_percentile_5 
-    if (all(c("mu_percentile_5","mu_percentile_95") %in% names(newdata@data))) 
-      newdata@data$interval90_mu <- newdata@data$mu_percentile_95 - newdata@data$mu_percentile_5 
-    if (all(c("prob_percentile_5","prob_percentile_95") %in% names(newdata@data))) 
-      newdata@data$interval90_prob <- newdata@data$prob_percentile_95 - newdata@data$prob_percentile_5 
-    if (all(c("Z_percentile_5","Z_percentile_95") %in% names(newdata@data))) 
-      newdata@data$interval90_Z <- newdata@data$Z_percentile_95 - newdata@data$Z_percentile_5 
-    
+
+    ## If the 5th and 95th percentiles are present, construct the 90% predictive interval width.
+    if (all(c("Y_percentile_5","Y_percentile_95") %in% names(newdata@data)))
+      newdata@data$interval90_Y <- newdata$Y_percentile_95 - newdata$Y_percentile_5
+    if (all(c("mu_percentile_5","mu_percentile_95") %in% names(newdata@data)))
+      newdata@data$interval90_mu <- newdata@data$mu_percentile_95 - newdata@data$mu_percentile_5
+    if (all(c("prob_percentile_5","prob_percentile_95") %in% names(newdata@data)))
+      newdata@data$interval90_prob <- newdata@data$prob_percentile_95 - newdata@data$prob_percentile_5
+    if (all(c("Z_percentile_5","Z_percentile_95") %in% names(newdata@data)))
+      newdata@data$interval90_Z <- newdata@data$Z_percentile_95 - newdata@data$Z_percentile_5
+
     ## Determine which quantities we are actually plotting
     QOI <- c("Y", "mu", "prob", "Z") # Possible Quantities of interest
     column_names <- paste0(rep(c("p_", "RMSPE_", "interval90_"), each = length(QOI)), QOI)
-    palette <- rep(c("Spectral", "BrBG", "BrBG"), each = length(QOI)) 
+    palette <- rep(c("Spectral", "BrBG", "BrBG"), each = length(QOI))
     keep <- column_names %in% names(newdata@data)
     column_names <- column_names[keep]
-    palette  <-  palette[keep] 
+    palette  <-  palette[keep]
   }
-  
-  ## plot() decides the palette for each graphic based on whether it is 
+
+  ## plot() decides the palette for each graphic based on whether it is
   ## the term is related to prediction or prediction uncertainty quantification.
-  ## This means that the palette argument cannot be provided by the user via 
-  ## ... argument. An exception is if the user has specified palette = "nasa", 
+  ## This means that the palette argument cannot be provided by the user via
+  ## ... argument. An exception is if the user has specified palette = "nasa",
   ## in which case we will use it for the prediction plots.
   user_args <- list(...)
   if (!is.null(user_args$palette) && user_args$palette == "nasa") {
     palette <- gsub(pattern = "Spectral", replacement = "nasa", palette)
   }
-  
-  ## Need to use do.call(), as it allows arguments to be provided as a list. 
+
+  ## Need to use do.call(), as it allows arguments to be provided as a list.
   ## Otherwise, there is no way to remove the original value of palette supplied
   ## by the user. Remove potential clashes:
   user_args$palette      <- NULL
   user_args$column_names <- NULL
-  args_list <- c(list(newdata = newdata, 
-                      column_names = column_names, 
-                      palette = palette), 
+  args_list <- c(list(newdata = newdata,
+                      column_names = column_names,
+                      palette = palette),
                  user_args)
   plots <- do.call(plot_spatial_or_ST, args_list)
-  
+
 
   ## First we determine if the user has predicted over the BAUs or over arbitrary
-  ## prediction regions. To assess this, we will test if newdata has the same 
+  ## prediction regions. To assess this, we will test if newdata has the same
   ## number of elements as the BAUs (this is not foolproof, but good enough).
   pred_over_BAUs <- length(object@BAUs) == length(newdata)
-  
+
   ## Edit labels
   split_column_names <- strsplit(column_names, "_")
   names(split_column_names) <- column_names
   for (i in column_names) {
     plots[[i]] <- plots[[i]] + .custom_lab(split_column_names[[i]], pred_over_BAUs)
   }
-  
-  ## Now plot the data. Note that here we use a simple call to plot_spatial_or_ST(), 
-  ## rather than via do.call(), because we have not constructed a custom palette 
+
+  ## Now plot the data. Note that here we use a simple call to plot_spatial_or_ST(),
+  ## rather than via do.call(), because we have not constructed a custom palette
   ## within this function.
   response_name <- all.vars(object@f)[1]
   if (is.data.frame(object@data)) {
@@ -343,18 +343,18 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
   } else if (is.list(object@data)) {
     if (any(sapply(object@data, function(d) is(d, "STIDF")))) {
       if (length(object@data) == 1) {
-        ## Plot the binned data over the BAUs. This is a bit difficult if any 
+        ## Plot the binned data over the BAUs. This is a bit difficult if any
         ## observations are associated with multiple BAUs; don't do it in this case.
         Cmat_dgT <- .as(object@Cmat, "dgTMatrix")
         if (!all(tabulate(Cmat_dgT@i + 1) == 1) || any(tabulate(Cmat_dgT@j + 1) > 1)) {
-         cat("The data set stored in object@data is of class STIDF. In this case, 
+         cat("The data set stored in object@data is of class STIDF. In this case,
               we normally use the binned data in object@Z to plot over the BAUs. However,
-              some observations are associated with multiple BAUs, or there 
-              are some BAUs associated with multiple observations (probably because 
+              some observations are associated with multiple BAUs, or there
+              are some BAUs associated with multiple observations (probably because
               average_in_BAU = FALSE). This complicates matters, so we will not plot the data.")
         } else {
-          Z <- binned_data(object) 
-          ## If NAs are present (i.e., some of the BAUs are unobserved), 
+          Z <- binned_data(object)
+          ## If NAs are present (i.e., some of the BAUs are unobserved),
           ## tell the user. This will illuminate a warning thrown later as well.
           if (any(is.na(Z))) {
             cat("To plot the STIDF data provided in the SRE object, we use the binned data in object@Z to plot over the BAUs. The unobserved BAUs (i.e., those that are not associated with any elements of object@Z) are given a value of NA.\n")
@@ -363,10 +363,10 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
           data_plots <- plot_spatial_or_ST(object@BAUs, response_name, ...)
         }
       } else {
-        ## We may have a combination of STIDF and STFDF, and we don't know which 
-        ## dataset each element of object@Z is associated with. This is a bit 
-        ## complicated, so we don't do it. 
-        cat("Multiple data sets were used in the analysis, and because at 
+        ## We may have a combination of STIDF and STFDF, and we don't know which
+        ## dataset each element of object@Z is associated with. This is a bit
+        ## complicated, so we don't do it.
+        cat("Multiple data sets were used in the analysis, and because at
                 least one is of class STIDF, we will not plot the data (we don't have a method for this).\n")
       }
     } else {
@@ -375,14 +375,14 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
   } else if (is(object@data, "SpatialPointsDataFrame") ||
              is(object@data, "SpatialPixelsDataFrame") ||
              is(object@data, "SpatialPolygonsDataFrame") ||
-             is(object@data, "STFDF")) { 
+             is(object@data, "STFDF")) {
     ## this shouldn't happen, but add just in case
     data_plots <- plot_spatial_or_ST(object@data, response_name, ...)
   } else {
     warning("Couldn't plot the data stored in object@data as the type was unrecognised.")
   }
-  
-  ## Give the data plots meaningful labels, and combine them with the prediction 
+
+  ## Give the data plots meaningful labels, and combine them with the prediction
   ## and uncertainty quantification plots.
   if (exists("data_plots")) {
     if (length(data_plots) == 1) {
@@ -394,19 +394,19 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
         data_plots[[i]] <- data_plots[[i]] + labs(colour = data_legend_name[i], fill = data_legend_name[i])
       }
     }
-    plots <- c(data_plots, plots) 
+    plots <- c(data_plots, plots)
   }
-  
+
   return(plots)
 }
 
 
 
 .custom_lab <- function(v, pred_over_BAUs) {
-    # v is a vector, with first entry indicating the type of plot we are making a 
-    ## label for, and the second entry indicating the quantity of interest    
+    # v is a vector, with first entry indicating the type of plot we are making a
+    ## label for, and the second entry indicating the quantity of interest
     type <- v[1]; x <- v[2]
-    
+
     ## Set the unicode character for the quantity of interest
     ## See https://en.wikipedia.org/wiki/List_of_Unicode_characters#Greek_and_Coptic
     # for the a list of unicode characters.
@@ -417,16 +417,16 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
     } else {
       process <-  bquote(bold(.(unicode))[P])
     }
-    
+
     expectation <- bquote(paste("E(", .(process), " | ", bold(Z), ", ", bold("\U03B8"), ")"))
-    
+
     ## Construct the labels
-    ## NB: Added white space to ensure no overlap between label and the fill box 
+    ## NB: Added white space to ensure no overlap between label and the fill box
     label <- if (type == "p") {
       # bquote(atop("Prediction    ", .(expectation))) # TODO add the following whitespace to expectation "    "
       top    <- "Prediction    "
       bottom <- bquote(paste(.(expectation) * "    "))
-      bquote(atop(.(top), .(bottom))) 
+      bquote(atop(.(top), .(bottom)))
     } else if (type == "RMSPE") {
         bquote(paste("RMSPE(" * .(expectation) * ", " * .(process),")        "))
     } else if (type == "interval90") {
@@ -434,7 +434,7 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
         bottom <- bquote(paste("interval width for " * .(process), "  "))
         bquote(atop(.(top), .(bottom)))
     }
-    
+
     return(labs(fill = label))
 }
 
@@ -444,68 +444,68 @@ setMethod("plot", signature(x = "SRE", y = "SpatialPolygonsDataFrame"), function
 
 #' @rdname plot_spatial_or_ST
 #' @export
-setMethod("plot_spatial_or_ST", signature(newdata = "STFDF"), 
-          function(newdata, column_names,  map_layer = NULL, 
-                   subset_time = NULL, palette = "Spectral", 
+setMethod("plot_spatial_or_ST", signature(newdata = "STFDF"),
+          function(newdata, column_names,  map_layer = NULL,
+                   subset_time = NULL, palette = "Spectral",
                    plot_over_world = FALSE, labels_from_coordnames = TRUE, ...) {
-      
-      ## Get the coordinate names 
+
+      ## Get the coordinate names
       coord_names <- coordnames(newdata)
-      
+
       ## Remove the coordinate corresponding to time. coord_names is just spatial.
       time_name <- coord_names[3]
       coord_names <- coord_names[1:2]
-      
-      if (!is.null(subset_time)) 
+
+      if (!is.null(subset_time))
           newdata <- newdata[, subset_time]
-                  
+
       plots <- .plot_spatial_or_ST_common(
-          newdata = newdata, column_names = column_names, coord_names = coord_names, 
-          time_name = time_name, map_layer = map_layer, 
-          palette = palette, plot_over_world = plot_over_world, 
+          newdata = newdata, column_names = column_names, coord_names = coord_names,
+          time_name = time_name, map_layer = map_layer,
+          palette = palette, plot_over_world = plot_over_world,
           labels_from_coordnames = labels_from_coordnames,
           ...
           )
       return(plots)
-})        
+})
 
 
 #' @rdname plot_spatial_or_ST
 #' @export
-setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPointsDataFrame"), 
-          function(newdata, column_names,  map_layer = NULL, 
-                   subset_time = NULL, palette = "Spectral", 
+setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPointsDataFrame"),
+          function(newdata, column_names,  map_layer = NULL,
+                   subset_time = NULL, palette = "Spectral",
                    plot_over_world = FALSE, labels_from_coordnames = TRUE, ...) {
-    
+
     ## Get the coordinate names, and time_name is NULL in the spatial setting
     coord_names <- coordnames(newdata)
     time_name <- NULL
-    
+
     .plot_spatial_or_ST_common(
-        newdata = newdata, column_names = column_names, coord_names = coord_names, 
-        time_name = time_name, map_layer = map_layer, 
-        palette = palette, plot_over_world = plot_over_world, 
-        labels_from_coordnames = labels_from_coordnames, 
+        newdata = newdata, column_names = column_names, coord_names = coord_names,
+        time_name = time_name, map_layer = map_layer,
+        palette = palette, plot_over_world = plot_over_world,
+        labels_from_coordnames = labels_from_coordnames,
         ...
     )
 })
 
 #' @rdname plot_spatial_or_ST
 #' @export
-setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPixelsDataFrame"), 
-          function(newdata, column_names,  map_layer = NULL, 
-                   subset_time = NULL, palette = "Spectral", 
+setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPixelsDataFrame"),
+          function(newdata, column_names,  map_layer = NULL,
+                   subset_time = NULL, palette = "Spectral",
                    plot_over_world = FALSE, labels_from_coordnames = TRUE,...) {
-            
+
             ## Get the coordinate names, and time_name is NULL in the spatial setting
             coord_names <- coordnames(newdata)
             time_name <- NULL
-            
+
             .plot_spatial_or_ST_common(
-              newdata = newdata, column_names = column_names, coord_names = coord_names, 
-              time_name = time_name, map_layer = map_layer, 
-              palette = palette, plot_over_world = plot_over_world, 
-              labels_from_coordnames = labels_from_coordnames, 
+              newdata = newdata, column_names = column_names, coord_names = coord_names,
+              time_name = time_name, map_layer = map_layer,
+              palette = palette, plot_over_world = plot_over_world,
+              labels_from_coordnames = labels_from_coordnames,
               ...
             )
           })
@@ -513,20 +513,20 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPixelsDataFrame"),
 
 #' @rdname plot_spatial_or_ST
 #' @export
-setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"), 
-          function(newdata, column_names,  map_layer = NULL, 
-                   subset_time = NULL, palette = "Spectral", 
+setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
+          function(newdata, column_names,  map_layer = NULL,
+                   subset_time = NULL, palette = "Spectral",
                    plot_over_world = FALSE, labels_from_coordnames = TRUE, ...) {
-            
+
             ## Get the coordinate names, and time_name is NULL in the spatial setting
             coord_names <- coordnames(newdata)
             time_name <- NULL
-            
+
             .plot_spatial_or_ST_common(
-              newdata = newdata, column_names = column_names, coord_names = coord_names, 
-              time_name = time_name, map_layer = map_layer, 
-              palette = palette, plot_over_world = plot_over_world, 
-              labels_from_coordnames = labels_from_coordnames, 
+              newdata = newdata, column_names = column_names, coord_names = coord_names,
+              time_name = time_name, map_layer = map_layer,
+              palette = palette, plot_over_world = plot_over_world,
+              labels_from_coordnames = labels_from_coordnames,
               ...
             )
           })
@@ -534,7 +534,7 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
 .plot_spatial_or_ST_common <- function(
     newdata, column_names, coord_names, time_name, map_layer, palette, plot_over_world, labels_from_coordnames, ...
 ) {
-  
+
     ## Inclusion of "-" characters can cause problems; convert to "_"
     column_names <- gsub("-", "_", column_names)
     names(newdata@data) <- gsub("-", "_", names(newdata@data))
@@ -543,33 +543,33 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
     tmp <- .classify_sp_and_convert_to_df(newdata)
     df      <- tmp$df
     sp_type <- tmp$sp_type
-    
-    ## Remove duplicate columns (can sometimes happen when we convert the Spatial* 
+
+    ## Remove duplicate columns (can sometimes happen when we convert the Spatial*
     ## object, e.g., if the coordinates are already present)
     df <- df[, unique(colnames(df))]
-    
+
     # ## Edit the time column so that the facets display t = ...
     # if(!is.null(time_name))
     # df[, time_name] <- factor(
     #    df[, time_name], ordered = TRUE,
     #    labels = paste(time_name, sort(unique(df[, time_name])), sep = " = ")
     # )
-    
-    if (length(palette) == 1) 
+
+    if (length(palette) == 1)
         palette <- rep(palette, length(column_names))
-  
-    
+
+
     ## Plot the requested columns
-    plots <- lapply(1:length(column_names), 
+    plots <- lapply(1:length(column_names),
                     function(i, x, y, ...) {
-                        .plot(df, x[i], coord_names, time_name, sp_type, map_layer, 
+                        .plot(df, x[i], coord_names, time_name, sp_type, map_layer,
                               y[i], plot_over_world, labels_from_coordnames, ...)
                     }, x = column_names, y = palette, ...)
-    
+
     suppressMessages( # suppress message about adding a new coordinate system
         if(plot_over_world) {
             plots <- lapply(plots, function(gg) {
-                (gg + .EmptyTheme_theme_only() + coord_map("mollweide")) %>%  
+                (gg + .EmptyTheme_theme_only() + coord_map("mollweide")) %>%
                     draw_world(inc_border = TRUE)
             })
         }
@@ -581,30 +581,30 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
 #### This function creates the individual plots.
 
 .plot <- function(df, column_name, coord_names, time_name, sp_type, map_layer, palette, plot_over_world, labels_from_coordnames, ...){
-    
+
     ## Basic plot
     if (!is.null(map_layer)) {
         if ("ggplot" %in% class(map_layer)) {
             ## Do it this way, because we cannot add map_layer to a ggplot object if
             ## map_layer is a gg object itself (e.g., if map_layer is a result from a
             ## call to ggmap(), which is often the case)
-            gg <- map_layer 
+            gg <- map_layer
             gg <- gg %+% df # change the default data to df
         } else {
-            ## If map_layer is not a ggplot object, it is probably a layer 
+            ## If map_layer is not a ggplot object, it is probably a layer
             ## (e.g., the result of a call to geom_polygon)
             gg <- ggplot(df) + map_layer
         }
-        
+
     } else {
-        gg <- ggplot(df) 
+        gg <- ggplot(df)
     }
     ## change/add the default aesthetics, and add some themes for nice plots
-    gg <- gg %+% aes_string(x = coord_names[1], y = coord_names[2]) + 
+    gg <- gg %+% aes_string(x = coord_names[1], y = coord_names[2]) +
         theme_bw() + coord_fixed()
-    
-    # ## If NAs are present, tell the user that we hard-code the colour/fill for 
-    # ## the pixels with NA values to be transparent 
+
+    # ## If NAs are present, tell the user that we hard-code the colour/fill for
+    # ## the pixels with NA values to be transparent
     # if (any(is.na(df[, column_name]))) {
     #   cat("NA values detected in the data, which will be transparent in the final plot. If you want them to show up in the plot, take the returned plot object and add a colour/fill scale with na.value equal to whatever colour you want.\n ")
     # }
@@ -615,29 +615,29 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
       colour_fn <- scale_colour_distiller(palette = palette, na.value = "transparent")
       fill_fn <- scale_fill_distiller(palette = palette, na.value = "transparent")
     }
-    
+
     ## Plot based on data type
     if (sp_type == "points") { # this is only for observation data
         gg <- gg + geom_point(aes_string(colour = column_name), ...) + colour_fn
     } else {
         if (sp_type == "pixels") {
-          ## geom_raster is typically faster, but it doesn't work if we are 
+          ## geom_raster is typically faster, but it doesn't work if we are
           ## plotting over the plane. For simplicity, just use geom_tile instead.
           if (plot_over_world) {
             gg <- gg + geom_tile(aes_string(fill = column_name), ...)
           } else {
             gg <- gg + geom_raster(aes_string(fill = column_name), ...)
           }
-          
+
         } else if (sp_type == "polygons") {
-            gg <- gg + geom_polygon(aes_string(group = "id", fill = column_name), ...) 
+            gg <- gg + geom_polygon(aes_string(group = "id", fill = column_name), ...)
         }
         gg <- gg + fill_fn
     }
-    
+
     if (!is.null(time_name))
         gg <- gg + facet_wrap(as.formula(paste("~", time_name)))
-    
+
     if (!labels_from_coordnames)
       gg <- gg + labs(x = expression(s[1]), y = expression(s[2]))
 
@@ -646,10 +646,10 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
 
 
 .classify_sp_and_convert_to_df <- function(newdata) {
-  
-    
+
+
     original_names <- names(newdata@data)
-    
+
     if (is(newdata, "SpatialPixelsDataFrame")) {
       df <- cbind(newdata@data, coordinates(newdata))
       sp_type <- "pixels"
@@ -665,15 +665,15 @@ setMethod("plot_spatial_or_ST", signature(newdata = "SpatialPolygonsDataFrame"),
         } else if (is(newdata@sp, "SpatialPixels")) {
             sp_type <- "pixels"
         }
-        df <- STFDF_to_df(newdata) 
+        df <- STFDF_to_df(newdata)
     } else {
         stop("Class of newdata is not recognised.")
     }
 
     ## check all original column names are present
     if (!all(original_names %in% names(df)))
-      warning("Some of the original names in newdata were not retained when newdata was coerced to a data.frame. 
+      warning("Some of the original names in newdata were not retained when newdata was coerced to a data.frame.
               This can sometimes happen if some of the column names contain '-', which get converted to '.'")
-    
+
     return(list(df = df, sp_type = sp_type))
 }
diff --git a/tests/testthat/test_basis.R b/tests/testthat/test_basis.R
index 6323888..aa9b60e 100644
--- a/tests/testthat/test_basis.R
+++ b/tests/testthat/test_basis.R
@@ -78,13 +78,18 @@ test_that("can average basis over polygons in plane", {
     HexPols_df <- SpatialPolygonsDataFrame(HexPols,
                                            cbind(over(HexPols,meuse.grid),
                                                  coordinates(HexPts)))
-    ## Generate basis functions
-    G <- auto_basis(manifold = plane(),data=meuse,nres = 2,prune=10,type = "Gaussian")
 
-        expect_true({eval_basis(G,coordinates(HexPts)); TRUE})
-        expect_true({eval_basis(G,HexPols_df); TRUE})
-        #plot(as.numeric(S1))
-        #lines(as.numeric(S2),col='red')
+    ## Generate regular basis functions with prune
+    G1 <- auto_basis(manifold = plane(),data=meuse,nres = 2,prune=10,type = "Gaussian")
+    expect_true({eval_basis(G1,coordinates(HexPts)); TRUE})
+    expect_true({eval_basis(G1,HexPols_df); TRUE})
+    #plot(as.numeric(S1))
+    #lines(as.numeric(S2),col='red')
+
+    ## Generate irregular basis functions, no prune, with fmesher
+    G2 <- auto_basis(manifold = plane(), data=meuse, nres = 2, type = "bisquare", regular = 0)
+    expect_true({eval_basis(G2,coordinates(HexPts)); TRUE})
+    expect_true({eval_basis(G2,HexPols_df); TRUE})
 })
 
 ## Deprecated: