#' fit_point_process
#'
#' @description Fit point process to randomize data
#'
#' @param pattern ppp object with point pattern
#' @param n_random Integer with number of randomizations.
#' @param process Character specifying which point process model to use.
#' Either \code{"poisson"} or \code{"cluster"}.
#' @param return_para Logical if fitted parameters should be returned.
#' @param return_input Logical if the original input data is returned.
#' @param simplify Logical if only pattern will be returned if \code{n_random = 1}
#' and \code{return_input = FALSE}.
#' @param verbose Logical if progress report is printed.
#'
#' @details
#' The functions randomizes the observed point pattern by fitting a point process to
#' the data and simulating \code{n_random} patterns using the fitted point process.
#' It is possible to choose between a Poisson process or a Thomas cluster process model.
#' For more information about the point process models, see e.g. Wiegand & Moloney (2014).
#'
#' @return rd_pat
#'
#' @examples
#' pattern_fitted <- fit_point_process(pattern = species_a, n_random = 39)
#'
#' @references
#' Plotkin, J.B., Potts, M.D., Leslie, N., Manokaran, N., LaFrankie, J.V.,
#' Ashton, P.S., 2000. Species-area curves, spatial aggregation, and habitat specialization
#' in tropical forests. Journal of Theoretical Biology 207, 81–99.
#' <https://doi.org/10.1006/jtbi.2000.2158>
#'
#' Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in
#' ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8
#'
#' @export
fit_point_process <- function(pattern, n_random = 1, process = "poisson", return_para = FALSE,
return_input = TRUE, simplify = FALSE, verbose = TRUE){
# check if n_random is >= 1
if (!n_random >= 1) {
stop("n_random must be >= 1.", call. = FALSE)
}
# unmark pattern
if (spatstat.geom::is.marked(pattern)) {
pattern <- spatstat.geom::unmark(pattern)
if (verbose) message("Unmarking provided input pattern.")
}
if (process == "poisson") {
result <- lapply(seq_len(n_random), function(x) {
simulated <- spatstat.random::runifpoint(n = pattern$n, win = pattern$window) # simulate poisson process
if (verbose) {
message("\r> Progress: n_random: ", x, "/", n_random, "\t\t", appendLF = FALSE)
}
return(simulated)
})
# calc parameters
if (return_para) {
number_points <- pattern$n
lambda <- spatstat.geom::intensity(pattern)
param_vec <- c(number_points = number_points,lambda = lambda)
}
} else if (process == "cluster") {
# fit cluster process
fitted_process <- spatstat.model::kppm(pattern, cluster = "Thomas",
statistic = "pcf",
statargs = list(divisor = "d",
correction = "best"),
method = "mincon", improve.type = "none")
# calc parameters
if (return_para) {
number_parents <- fitted_process$clustpar[["kappa"]] * spatstat.geom::area(pattern$window)
number_points <- fitted_process$mu
cluster_area <- fitted_process$clustpar[["scale"]] ^ 2 * pi
param_vec <- c(number_parents = number_parents, number_points = number_points,
cluster_area = cluster_area)
}
result <- lapply(seq_len(n_random), function(x) {
# simulate clustered pattern
simulated <- spatstat.model::simulate.kppm(fitted_process, window = pattern$window,
nsim = 1, drop = TRUE)
# remove points because more points in simulated
if (pattern$n < simulated$n) {
# difference between patterns
difference <- simulated$n - pattern$n
# id of points to remove
remove_points <- sample(x = seq_len(simulated$n), size = difference)
# remove points
simulated <- simulated[-remove_points]
# add points because less points in simulated
} else if (pattern$n > simulated$n) {
# difference between patterns
difference <- pattern$n - simulated$n
# create missing points
missing_points <- spatstat.random::runifpoint(n = difference, win = pattern$window,
nsim = 1, drop = TRUE)
# add missing points to simulated
simulated <- spatstat.geom::superimpose(simulated, missing_points,
W = pattern$window)
}
if (verbose) {
message("\r> Progress: n_random: ", x, "/", n_random, "\t\t", appendLF = FALSE)
}
return(simulated)
})
} else {
stop("Please select either 'poisson' or 'cluster'.", call. = FALSE)
}
# set param to NA
if (!return_para) param_vec <- NA
# set names
names(result) <- paste0("randomized_", seq_len(n_random))
# combine to one list
result <- list(randomized = result, observed = pattern, method = "fit_point_process()",
energy_df = NA, stop_criterion = NA, iterations = NA, param = param_vec)
# add param
if (return_para) {
result$param <- param_vec
}
# set class of result
class(result) <- "rd_pat"
# remove input if return_input = FALSE
if (!return_input) {
# set observed to NA
result$observed <- NA
# check if output should be simplified
if (simplify) {
# not possible if more than one pattern is present
if (n_random > 1) {
warning("'simplify = TRUE' not possible for 'n_random > 1'.",
call. = FALSE)
# only one random pattern is present that should be returend
} else if (n_random == 1) {
result <- result$randomized[[1]]
}
}
# return input if return_input = TRUE
} else {
# return warning if simply = TRUE because not possible if return_input = TRUE (only verbose = TRUE)
if (simplify) {
warning("'simplify = TRUE' not possible for 'return_input = TRUE'.", call. = FALSE)
}
}
# write result in new line if progress was printed
if (verbose) {
message("\r")
}
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.