Nothing
#' @title Spatial kernel density estimate
#' @description A weighted or unweighted Gaussian Kernel Density estimate
#' for point spatial data
#'
#' @param x sf POINT object
#' @param y Optional values, associated with x coordinates,
#' to be used as weights
#' @param bw Distance bandwidth of Gaussian Kernel, must be units
#' of projection
#' @param ref A terra SpatRaster, sf class object or c[xmin,xmax,ymin,ymax]
#' vector to estimate the kde extent
#' @param res Resolution of raster when ref not SpatRaster
#' @param standardize Standardize results to 0-1 (FALSE/TRUE)
#' @param scale.factor Numeric scaling factor for the KDE (defaults to 10000),
#' to account for very small estimate values
#' @param mask (TRUE/FALSE) mask resulting raster if ref is provided
#' as a SpatRaster
#'
#' @details
#' Please note that ks methods for estimation has been reverted to the Gussian method proposed
#' in Venables & Ripley (2002). There was not enought evendence that the Chacon & Duong (2018)
#' multivariate method(s) for bandwidth selection and kernal estimation were suitable for
#' spatial random fields.
#'
#' @return a terra SpatRaster class object containing kernel density estimate
#'
#' @author Jeffrey S. Evans <jeffrey_evans@@tnc.org>
#'
#' @references
#' Duong, T. & Hazelton, M.L. (2005) Cross-validation bandwidth matrices for multivariate
#' kernel density estimation. Scandinavian Journal of Statistics, 32, 485-506.
#'
#' Wand, M.P. & Jones, M.C. (1994) Multivariate plug-in bandwidth selection. Computational
#' Statistics, 9, 97-116.
#'
#' Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S.
#' Fourth edition. Springer.
#'
#' @examples
#' \donttest{
#' library(sf)
#' library(terra)
#'
#' data(meuse, package = "sp")
#' meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992,
#' agr = "constant")
#'
#' # Unweighted KDE (spatial locations only) with 40m resoultion
#' pt.kde <- sf.kde(x = meuse, bw = 1000, standardize = TRUE, res=40)
#' plot(pt.kde, main="Unweighted kde")
#' plot(st_geometry(meuse), pch=20, col="red", add=TRUE)
#'
#' # cadmium weighted KDE usign extent with 40m resoultion and 500m and 1000m bw
#' cadmium.kde.500 <- sf.kde(x = meuse, y = meuse$cadmium, res=40,
#' bw = 500, standardize = TRUE)
#' cadmium.kde.1000 <- sf.kde(x = meuse, y = meuse$cadmium, res=40,
#' bw = 1000, standardize = TRUE)
#' plot(c(cadmium.kde.500, cadmium.kde.1000))
#'
#' }
#'
#' @export sp.kde sf.kde
#' @aliases sp.kde
sf.kde <- function(x, y = NULL, bw = NULL, ref = NULL, res = NULL,
standardize = FALSE, scale.factor = 10000,
mask = FALSE) {
if(missing(x))
stop("x argument must be provided")
if(!inherits(x, c("sf", "sfc") ))
stop(deparse(substitute(x)), " must be a sf, or sfc object")
if(unique(as.character(sf::st_geometry_type(x))) != "POINT")
stop(deparse(substitute(x)), " must be single-part POINT geometry")
if(is.null(scale.factor)) scale.factor = 1
ref.flag = inherits(ref, "SpatRaster")
if(is.null(res) && ref.flag == TRUE) {
res <- terra::res(ref)
}
if(is.null(res) && ref.flag == FALSE) {
res <- 100
message("resoultion not defined, defaulting to 100")
}
if(is.null(ref)) {
ref <- terra::rast(terra::ext(x), resolution = res,
crs=terra::crs(x))
} else if(inherits(ref, "numeric")) {
if(length(ref) != 4)
stop("Need xmin, xmax, ymin, ymax bounding coordinates")
ref <- terra::rast(terra::ext(ref), resolution = res,
crs=terra::crs(x))
} else {
if(!inherits(ref, "SpatRaster"))
stop(deparse(substitute(ref)), " must be a terra SpatRast object")
if(terra::res(ref)[1] != res)
message("reference raster defined, res argument is being ignored")
}
# weighted kde function, modification of MASS::kde2d
fhat <- function (x, y, h, w, n = 25, lims = c(range(x), range(y))) {
nx <- length(x)
if (length(y) != nx)
stop("data vectors must be the same length")
if (length(w) != nx & length(w) != 1)
stop("weight vectors must be 1 or length of data")
if (missing(h)) {
h <- c(MASS::bandwidth.nrd(x), MASS::bandwidth.nrd(y))
} else {
h <- rep(h, length.out = 2L)
}
if (any(h <= 0)) stop("bandwidths must be strictly positive")
if (missing(w)) { w <- numeric(nx) + 1 }
gx <- seq(lims[1], lims[2], length = n[1])
gy <- seq(lims[3], lims[4], length = n[2])
h <- h/4
ax <- outer(gx, x, "-") / h[1]
ay <- outer(gy, y, "-") / h[2]
z <- ( matrix(rep(w, n[1]), nrow = n[1], ncol = nx, byrow = TRUE) *
matrix(stats::dnorm(ax), n[1], nx) ) %*% t(matrix(stats::dnorm(ay), n[2], nx)) /
( sum(w) * h[1] * h[2] )
return(list(x = gx, y = gy, z = z))
}
if(is.null(bw)){
bwf <- function(x){
r <- stats::quantile(x, c(0.25, 0.75))
h <- (r[2] - r[1])/1.34
4 * 1.06 * min(sqrt(stats::var(x)), h) * length(x)^(-1/5)
}
bw <- c(bwf(sf::st_coordinates(x)[,1]),
bwf(sf::st_coordinates(x)[,2]))
message("Using ", round(bw[1],3), ", ", round(bw[2],3), " for bandwidth", "\n")
} else {
bw <- c(bw,bw)
}
n <- c(terra::nrow(ref), terra::ncol(ref))
if(!is.null(y)) {
message("\n","calculating weighted kde","\n")
k <- fhat(sf::st_coordinates(x)[,1], sf::st_coordinates(x)[,2], w = y,
h = bw, n = n, lims = as.vector(terra::ext(ref)) )
} else {
message("\n","calculating unweighted kde","\n")
k <- MASS::kde2d(sf::st_coordinates(x)[,1], sf::st_coordinates(x)[,2], h = bw,
n = n, lims = as.vector(terra::ext(ref)) )
}
k$z <- k$z * scale.factor
if( standardize == TRUE ) { k$z <- (k$z - min(k$z)) / (max(k$z) - min(k$z)) }
pts <- data.frame(expand.grid(x=k$x, y=k$y),
z=round(as.vector(array(k$z,length(k$z))) *
scale.factor, 10))
kde.est <- terra::rast(pts, type="xyz", extent = terra::ext(ref) )
if(mask == TRUE) {
kde.est <- terra::mask(kde.est, ref)
}
terra::crs(kde.est) <- terra::crs(x)
return( kde.est )
}
sp.kde = sf.kde
# sf.kde <- function(x, y = NULL, bw = NULL, ref = NULL, res = NULL,
# standardize = FALSE, scale.factor = NULL,
# mask = FALSE) {
# if(missing(x))
# stop("x argument must be provided")
# ref.flag = inherits(ref, "SpatRaster")
# if(!inherits(x, c("sf", "sfc") ))
# stop(deparse(substitute(x)), " must be a sf, or sfc object")
# if(unique(as.character(sf::st_geometry_type(x))) != "POINT")
# stop(deparse(substitute(x)), " must be single-part POINT geometry")
# if(is.null(ref)) {
# if(!is.null(res)) {
# ref <- terra::rast(terra::ext(x), resolution = res)
# } else {
# ref <- terra::rast(terra::ext(x))
# message("defaulting to ", terra::res(ref)[1], "x", terra::res(ref)[2], " cell resolution")
# }
# }
# if(inherits(ref, "numeric")) {
# if(length(ref) != 4)
# stop("Need xmin, xmax, ymin, ymax bounding coordinates")
# if(!is.null(res)) {
# ref <- terra::rast(terra::ext(ref), resolution = res)
# } else {
# ref <- terra::rast(terra::ext(ref))
# message("defaulting to ", terra::res(ref)[1], "x", terra::res(ref)[2], " cell resolution")
# }
# } else {
# if(!inherits(ref, "SpatRaster"))
# stop(deparse(substitute(ref)), " must be a terra SpatRast object")
# }
# n <- c(terra::nrow(ref), terra::ncol(ref))
# if(!is.null(bw)) {
# message("Using specified bandwidth: "); print(bw)
# } else {
# if(is.null(y)) {
# bw = suppressWarnings(ks::hpi(sf::st_coordinates(x)[,1:2]))
# message("Unweighted automatic bandwidth: "); print(bw)
# } else {
# bw = suppressWarnings(ks::Hscv.diag(cbind(sf::st_coordinates(x)[,1:2],y)))
# message("Weighted automatic CV bandwidth: "); print(bw)
# }
# }
# if(!is.null(y)) {
# message("\n","calculating weighted kde","\n")
# kde.est <- suppressWarnings(
# terra::rast(matrix(ks::kde(sf::st_coordinates(x)[,1:2],
# h=bw, eval.points=terra::xyFromCell(ref, 1:terra::ncell(ref)),
# gridsize=n, w = y, density=TRUE)$estimate,
# nrow=n[1], ncol=n[2], byrow=TRUE),
# extent=terra::ext(ref)) )
# } else {
# message("\n","calculating unweighted kde","\n")
# kde.est <- suppressWarnings(
# terra::rast(matrix(ks::kde(sf::st_coordinates(x)[,1:2],
# h=bw, eval.points=terra::xyFromCell(ref, 1:terra::ncell(ref)),
# gridsize=n, density=TRUE)$estimate,
# nrow=n[1], ncol=n[2], byrow=TRUE),
# extent=terra::ext(ref)) )
# }
# if(!is.null(scale.factor)) kde.est <- kde.est * scale.factor
# if( standardize == TRUE ) { kde.est <- kde.est /
# terra::global(kde.est, "max", na.rm=TRUE)[,1] }
# if(mask) {
# if(!ref.flag) {
# message("Since a raster was not used as ref, there is nothing to mask")
# } else {
# kde.est <- terra::mask(kde.est, ref)
# }
# }
# terra::crs(kde.est) <- terra::crs(x)
# return( kde.est )
# }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.