diff --git a/NAMESPACE b/NAMESPACE
index 385af71..0f4b2d0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,7 @@
 
 S3method(as.data.frame,Basis)
 S3method(as.data.frame,TensorP_Basis)
+export(BAUs_from_points)
 export(Basis_as_data.frame)
 export(EmptyTheme)
 export(Euclid_dist)
diff --git a/R/AllGeneric.R b/R/AllGeneric.R
index c7e609d..c05e8d7 100644
--- a/R/AllGeneric.R
+++ b/R/AllGeneric.R
@@ -171,7 +171,24 @@ TensorP_Basis_as_data.frame <- setAs("TensorP_Basis", "data.frame",
 #' @rdname Basis_data.frame
 setGeneric("data.frame<-", function(x, value) standardGeneric("data.frame<-"))
 
-#### NOT EXPORTED ####
+
+#' @title Creates pixels around points
+#' @description Takes a SpatialPointsDataFrame and converts it into SpatialPolygonsDataFrame by constructing a tiny (within machine tolerance) BAU around each SpatialPoint.
+#' @param sp_obj object of class \code{SpatialPointsDataFrame}
+#' @details This function is there to allow users to mimic standard geospatial analysis where BAUs are not used. Since \code{FRK} is build on the concept of a BAU, this function constructs tiny BAUs around the observation and prediction locations which can be subsequently passed on to the functions \code{SRE} and \code{FRK}. With \code{BAUs_from_points}, the user supplies both the data and prediction locations accompanied with covariates.
+#' @export
+#' @examples
+#' library(sp)
+#' opts_FRK$set("parallel",0L)
+#' df <- data.frame(x = rnorm(10),
+#'                  y = rnorm(10))
+#' coordinates(df) <- ~x+y
+#' BAUs <- BAUs_from_points(df)
+setGeneric("BAUs_from_points", function(obj) standardGeneric("BAUs_from_points"))
+
+########################
+#### NOT EXPORTED ######
+########################
 
 #' @title Automatic BAU generation
 #' @noRd
diff --git a/R/FRK_wrapper.R b/R/FRK_wrapper.R
index be3c1a1..7dad867 100644
--- a/R/FRK_wrapper.R
+++ b/R/FRK_wrapper.R
@@ -167,7 +167,7 @@ FRK <- function(f,                     # formula (compulsory)
              please supply a field 'std' in the data objects.")
     if(!(is.null(BAUs))) {
         if(!(is(BAUs,"SpatialPolygonsDataFrame") | is(BAUs,"SpatialPixelsDataFrame") | is(BAUs,"STFDF")))
-            stop("BAUs should be a SpatialPolygonsDataFrame or a STFDF object")
+            stop("BAUs should be a SpatialPolygonsDataFrame, SpatialPixelsDataFrame, or a STFDF object")
         if(!all(sapply(data,function(x) identical(proj4string(x), proj4string(BAUs)))))
             stop("Please ensure all data items and BAUs have the same coordinate reference system")
         if(!(all(BAUs$fs >= 0)))
diff --git a/R/geometryfns.R b/R/geometryfns.R
index f4d90e5..eb2a8e5 100644
--- a/R/geometryfns.R
+++ b/R/geometryfns.R
@@ -982,6 +982,51 @@ SpatialPolygonsDataFrame_to_df <- function(sp_polys,vars = names(sp_polys)) {
     df_polys
 }
 
+
+## Create very small square polygon BAUs around SpatialPoints
+setMethod("BAUs_from_points",signature(obj = "SpatialPoints"),
+          function(obj) {
+
+    sp_obj_pols <- NULL                 # Initialise polygons
+    offset <- 10*.Machine$double.eps    # Side size of BAU (very small)
+    cnames <- coordnames(obj)           # coordinate names
+    coords <- coordinates(obj)          # coordinates of SpatialPoints
+
+
+    ## Generate the Bottom Left, Bottom Right, Top Right, and Top Left, corners of the BAUs
+    BL <- data.frame(X1 = coords[,1] - offset, X2 = coords[,2] - offset, id = 1:length(obj))
+    BR <- data.frame(X1 = coords[,1] + offset, X2 = coords[,2] - offset, id = 1:length(obj))
+    TR <- data.frame(X1 = coords[,1] + offset, X2 = coords[,2] + offset, id = 1:length(obj))
+    TL <- data.frame(X1 = coords[,1] - offset, X2 = coords[,2] + offset, id = 1:length(obj))
+
+    ## Interleave them appropriate so they form polygon paths and set names
+    sp_obj_pols <- .interleave(BL,BR,TR,TL)
+    names(sp_obj_pols) <- c(cnames,"id")
+
+    ## Now create polygons from the above paths, and keep same projection
+    sp_obj_pols <- df_to_SpatialPolygons(sp_obj_pols,coords=cnames,keys="id",
+                                         proj = CRS(proj4string(obj)))
+
+    ## We assign the centroid of the BAU to the data object
+    df_data <- as.data.frame(coords)
+
+    ## If data points had other variables, add them aswell
+    if(is(obj,"SpatialPointsDataFrame"))
+        df_data <- cbind(df_data,obj@data)
+
+    ## Ensure the row names are the same and construct the SpatialPolygonsDataFrame BAUs
+    row.names(df_data) <- row.names(sp_obj_pols)
+    sp_obj_pols <- SpatialPolygonsDataFrame(sp_obj_pols,data = df_data)
+})
+
+## Create very small square polygon BAUs around SpatialPoints
+setMethod("BAUs_from_points",signature(obj = "ST"),
+          function(obj) {
+             print("BAUs from points for space-time data not yet implemented. Please contact
+                   the package maintainer.")
+          })
+
+
 ###########################################
 ########## Not exported ###################
 ###########################################
@@ -1690,3 +1735,25 @@ process_isea3h <- function(isea3h,resl) {
    } else { LHS <- dep_var }                                 # otherwise it's just the dep var
    newf <- formula(paste0(LHS,"~1"))                         # now return formula without covariates
 }
+
+## Interleave data frames
+.interleave <- function(...) {
+
+    ## Extract data frames
+    dfs <- list(...)
+
+    ## Basic checks -- all data frames, same names and all same length
+    stopifnot(all(sapply(dfs,function(x) is(x,"data.frame"))))
+    stopifnot(all(sapply(dfs,function(x) names(x) == names(dfs[[1]]))))
+    stopifnot(all(sapply(dfs,function(x) nrow(x) == nrow(dfs[[1]]))))
+
+    ndfs <- length(dfs)   # number of data frames
+    n <- nrow(dfs[[1]])   # number of rows
+
+    stacked <- do.call("rbind",dfs)                       # stack all together
+    idx <- rep(1:n, each = ndfs) + (0:(ndfs-1)) * n       # interleaving indices
+    interleaved <- stacked[idx,]                          # re-order appropriately
+    row.names(interleaved) <- 1:nrow(interleaved)         # reset row names
+    interleaved                                           # return
+
+}
diff --git a/man/BAUs_from_points.Rd b/man/BAUs_from_points.Rd
new file mode 100644
index 0000000..d044fc5
--- /dev/null
+++ b/man/BAUs_from_points.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/AllGeneric.R
+\name{BAUs_from_points}
+\alias{BAUs_from_points}
+\title{Creates pixels around points}
+\usage{
+BAUs_from_points(obj)
+}
+\arguments{
+\item{sp_obj}{object of class \code{SpatialPointsDataFrame}}
+}
+\description{
+Takes a SpatialPointsDataFrame and converts it into SpatialPolygonsDataFrame by constructing a tiny (within machine tolerance) BAU around each SpatialPoint.
+}
+\details{
+This function is there to allow users to mimic standard geospatial analysis where BAUs are not used. Since \code{FRK} is build on the concept of a BAU, this function constructs tiny BAUs around the observation and prediction locations which can be subsequently passed on to the functions \code{SRE} and \code{FRK}. With \code{BAUs_from_points}, the user supplies both the data and prediction locations accompanied with covariates.
+}
+\examples{
+library(sp)
+opts_FRK$set("parallel",0L)
+df <- data.frame(x = rnorm(10),
+                 y = rnorm(10))
+coordinates(df) <- ~x+y
+BAUs <- BAUs_from_points(df)
+}
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
new file mode 100644
index 0000000..cbe0da2
Binary files /dev/null and b/tests/testthat/Rplots.pdf differ
diff --git a/tests/testthat/test_BAUs.R b/tests/testthat/test_BAUs.R
index e76acbb..8c6bf13 100644
--- a/tests/testthat/test_BAUs.R
+++ b/tests/testthat/test_BAUs.R
@@ -226,3 +226,19 @@ test_that("SpaceTime_BAUs",{
     expect_equal(attr(space_time_grid@time,"tzone"),attr(STobj1@time,"tzone"))
 })
 
+
+test_that("Point from BAUs works",{
+    library(sp)
+    dat <- data.frame(x = rnorm(100),
+                      y = rnorm(100))
+    coordinates(dat) <- ~x+y
+    BAUs <- BAUs_from_points(dat)
+
+    expect_is(BAUs,"SpatialPolygonsDataFrame")
+    expect_equal(length(BAUs),100)
+
+    dat$z <- rnorm(100)
+    BAUs <- BAUs_from_points(dat)
+    expect_is(BAUs,"SpatialPolygonsDataFrame")
+    expect_equal(length(BAUs),100)
+})