Nothing
#' Kernel Density Estimation
#'
#' KDE for spatial data. The algorithm is heavily inspired by
#' \href{https://github.com/qgis/QGIS/blob/b3d2619976a69d7fb67b884492da491dfaba287c/src/analysis/raster/qgskde.cpp}{Heatmap tool}
#' in QGIS. The help for QGIS tools is provided \href{https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/qgis/interpolation.html#heatmap-kernel-density-estimation}{at the QGIS website}.
#' The a tutorial is provided \href{https://grindgis.com/software/heat-map-using-qgis}{here}.
#'
#' @details
#' \code{grid} parameter specifies output of the function. KDE is calculated on the specified \code{grid}.
#' If grid is \code{\link[raster]{Raster-class}} then outcome is also \code{\link[raster]{Raster-class}}.
#' If grid is \code{\link[sf]{sf}} \code{data.frame} then outcome is also \code{\link[sf]{sf}} \code{data.frame}.
#'
#'
#' @param points \code{\link[sf]{sf}} \code{data.frame} containing only POINTS.
#' @param band_width \code{numeric} specifying the band width for KDE.
#' @param decay \code{numeric} specifying the decay parameter for \code{"triangular"} kernel. For
#' other kernels besides \code{"triangular"} the parameter is not used.
#' @param kernel \code{character} specifying type of kernel to use. Available implemented kernels are
#' \code{"uniform", "quartic", "triweight", "epanechnikov", "triangular"}. Default is \code{"quartic"} and if
#' unknown kernel name is used it falls back to the default value.
#' @param scaled \code{logical} specifying if the output values should be scaled. Default value is
#' \code{FALSE}.
#' @param weights \code{numeric} vector of weights for individual \code{points}.
#' @param grid either \code{\link[sf]{sf}} \code{data.frame} (outcome of function
#' \code{\link{create_grid_rectangular}} or \code{\link{create_grid_hexagonal}}) or
#' \code{\link[raster]{Raster-class}} (outcome of function \code{\link{create_raster}}).
#' Does not have to be specified if \code{cell_size} is set.
#' @param cell_size \code{numeric} specifying the distance for equal spaced points. Must be
#' higher than 0. Can be left out if \code{grid} is provided as \code{grid} is used instead.
#' The code used to generate grid is \code{\link{create_grid_rectangular}(points, cell_size, band_width)}.
#' @param quiet Should printing of progress bar be suppressed? Default `FALSE`.
#'
#' @return either \code{\link[sf]{sf}} \code{data.frame} or \code{\link[raster]{Raster-class}}
#' depending on class of \code{grid} parameter.
#' @export
#'
#' @importFrom sf st_is_longlat st_bbox st_geometry st_union st_convex_hull st_buffer st_make_grid
#' @importFrom sf st_sf st_coordinates st_geometry_type st_centroid
#' @importFrom dplyr mutate
#' @importFrom glue glue glue_collapse
#' @importFrom rlang .data is_double is_integer
#' @importFrom vctrs vec_cast
#'
#' @useDynLib SpatialKDE, .registration = TRUE
#'
#' @examples
#' library(sf)
#' nc <- st_read(system.file("shape/nc.shp", package = "sf")) %>% st_transform(32031)
#' grid <- create_grid_hexagonal(nc, cell_size = 100000)
#' points <- st_sample(nc, 500) %>% st_as_sf()
#' kde_estimate_grid <- kde(points, band_width = 150000, grid = grid)
#' raster <- create_raster(nc, cell_size = 100000)
#' kde_estimate_raster <- kde(points, band_width = 150000, grid = raster)
#'
kde <- function(points,
band_width,
decay = 1,
kernel = c("quartic", "uniform", "triweight", "epanechnikov", "triangular"),
scaled = FALSE,
weights = c(),
grid,
cell_size,
quiet = FALSE) {
kernel <- match.arg(kernel)
.validate_sf(points)
.validate_projected(points)
.validate_points(points)
if (length(weights) == 0) {
weights <- rep(1, nrow(points))
}
if (is_integer(weights, n = nrow(points))) {
weights <- vctrs::vec_cast(weights, double())
}
if (!is_double(weights, n = nrow(points), finite = TRUE)) {
stop(glue::glue(
"All values of `weights` must be numerical and finite vector (no `NA`s, `Inf` or `-Inf`).",
"The length of the vector must be equal to number of rows in `points`.",
"Length weights is `{length(weights)}` and number of rows in points is `{nrow(points)}`."
))
}
if (missing(grid) & missing(cell_size)) {
stop("Both variables `grid` and `cellsize` are not specified. Don't know how to create grid for KDE estimation.")
}
if (!missing(cell_size)) {
.validate_cellsize(cell_size)
}
if (missing(grid)) {
grid <- create_grid_rectangular(points, cell_size, band_width)
}
.validate_bandwidth(band_width)
kde_calculated <- .kde(
grid = grid,
points = points,
band_width = band_width,
decay = decay,
kernel = kernel,
scaled = scaled,
weights = weights,
quiet = quiet
)
return(kde_calculated)
}
.kde <- function(grid,
points,
band_width,
decay,
kernel,
scaled,
weights,
quiet = FALSE) {
UseMethod(".kde")
}
#' @importFrom sf st_centroid st_coordinates st_crs st_transform
#' @importFrom dplyr mutate
#' @importFrom glue glue
.kde.sf <- function(grid,
points,
band_width,
decay,
kernel,
scaled,
weights,
quiet) {
.validate_sf(grid)
.validate_projected(grid)
cells_number <- nrow(grid)
.warn_long_calculation(cells_number)
if (all(unique(st_geometry_type(grid)) == "POINT")) {
grid_points <- grid
} else {
if (check_message_print()) {
message("Using centroids instead of provided `grid` geometries to calculate KDE estimates.")
}
suppressWarnings(
grid_points <- grid %>%
sf::st_centroid(of_largest_polygon = TRUE)
)
}
if (sf::st_crs(points) != sf::st_crs(grid_points)) {
message(glue::glue(
"`points` are transformed into `grid` CRS to make coordinates match. ",
"You may want to do this manually using `st_transform()` function to ",
"have better control over the process."
))
points <- points %>%
sf::st_transform(crs = st_crs(grid_points))
}
kde_values <- kde_estimate(sf::st_coordinates(grid_points),
sf::st_coordinates(points),
bw = band_width,
kernel = kernel,
scaled = scaled,
decay = decay,
weights = weights,
showProgressBar = !quiet
)
grid <- grid %>%
dplyr::mutate(kde_value = kde_values)
grid
}
#' @importFrom raster xyFromCell values
#' @importFrom sf st_coordinates
#' @importFrom methods setMethod
setMethod(
".kde",
"RasterLayer",
function(grid,
points,
band_width,
decay,
kernel,
scaled,
weights,
quiet) {
.validate_raster_projected(grid)
cells_number <- length(grid)
.warn_long_calculation(cells_number)
kde_values <- kde_estimate(raster::xyFromCell(grid, 1:cells_number),
sf::st_coordinates(points),
bw = band_width,
kernel = kernel,
scaled = scaled,
decay = decay,
weights = weights,
showProgressBar = !quiet
)
raster::values(grid) <- kde_values
grid
}
)
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.