R/pcf_calc.R
In jnghiem/bfasttools: Tools for BFAST spatial analysis
Defines functions
pcf_calc
Documented in
pcf_calc
#' Compute pair correlation function for breakpoint distribution
#'
#' This function computes the pair correlation function (PCF) for the point
#' pattern of breakpoints. The PCF estimates \emph{g} for a given radius
#' distance \emph{r}, where \emph{g} represents the degree of clustering.
#'
#' This function is a convenient wrapper for \code{spatstat::pcf()} and
#' \code{spatstat::envelope()}. Most default values in \code{pcf} and
#' \code{envelope} are accepted. However, the translation correction is used
#' here.
#'
#' If \code{level} is not set to \code{NA}, a simulation envelope will be
#' computed according to the number of simulations needed to achieve
#' \code{level}. However, the simulation process is very time-consuming so
#' setting \code{level=NA} may be best choice for reasonable processing times.
#'
#' The returned data frame has the following fields: \enumerate{ \item
#' \code{year}: The year of the observation. \item \code{month}: The month of
#' the observation. \item \code{type}: The type of statistic for \emph{g}. \item
#' \code{g}: The value of \emph{g} for the given \code{type}. \item \code{dist}:
#' The radius distance. \item \code{upper.bound, lower.bound}: The upper/lower
#' bound of the simulation envelope at distance \code{dist}. \item \code{nsim}:
#' The number of simulations to produce the envelope. \item \code{sig.level}:
#' The significance level of the pointwise Monte Carlo test. }
#'
#' @param bfast_in A data frame generated by \code{bfastSpatial()}.
#' @param level A number indicating the significance level for the pointwise
#' Monte Carlo test. See the documentation for \code{spatstat::envelope()} for
#' more details. If \code{level=NA}, no envelope will be simulated.
#' @param boundary A SpatialPolygons object defining the boundary of the PCF
#' analysis. If this argument is not supplied, all records in \code{bfast_in}
#' will be used.
#' @param template A Raster object with the same resolution, extent, and
#' projection as the raster dataset used to generate \code{bfast_in}.
#' Alternatively, an XML file generated by \code{create_raster_metadata()}.
#' @param resample_dist An optional number indicating the square dimension to
#' which the breakpoint pattern raster should be resampled. This argument is
#' in the same units as the projection of \code{template}. If not supplied, no
#' resampling will occur and points will be generated based on the resolution
#' of \code{template}.
#' @param max.dist The maximum radius distance to which the PCF calculation
#' should be performed. If not supplied, the default is the diagonal distance
#' across \code{boundary}.
#' @param output_directory An optional character file path to a directory to
#' which the PCF calculations for each year and month should be written in
#' \code{.rds} format. The radius distances for which the PCF was computed
#' will also be written here.
#' @param mc.cores A numeric indicating the number of cores to be used in
#' parallel computing.
#' @return A list of length 2. The first item is a data frame summarizing the
#' PCF results. The second item is the vector of radius distances.
#' @examples
#' \dontrun{
#' pcf_calc(bf_df, NA, rgdal::readOGR(dsn="C:/Desktop/shapefiles", layer="Mojave"), template, 2000, "C:/Desktop/PCF_files", 6)
#' }
#' @importFrom spatstat owin Window pcf envelope
#' @importFrom foreach %dopar% foreach
#' @importFrom maptools as.ppp.SpatialPointsDataFrame
#' @export
pcf_calc <- function(bfast_in, level=0.05, boundary, template, resample_dist=NULL, max.dist=NULL, output_directory=NULL, mc.cores=1) {
template <- create_template_raster(template)
if (!is.null(boundary)) {
cells <- extract(template, boundary, cellnumbers=TRUE, df=TRUE)[,"cell"]
} else {
cells <- 1:ncell(template)
}
results_cleaned <- bfast_in %>%
clean_bfast() %>%
filter(no_cell %in% cells, !is.na(start_date))
ext <- extent(boundary)
window <- owin(xrange=ext[1:2], yrange=ext[3:4])
brks_only <- filter(results_cleaned, brk==1)
years <- results_cleaned[,"year"]; yr.range <- min(years):max(years)
simulations <- round((2/level)-1)
if (is.null(max.dist)) {
max.dist <- ceiling(sqrt(diff(ext[1:2])^2+diff(ext[3:4])^2)/2)
}
if (!is.null(resample_dist)) {
resample_ras <- raster(ext, resolution=resample_dist, crs=proj4string(template))
radii <- seq(0, max.dist, by=signif(res(resample_ras)[1], 2))
} else {
radii <- seq(0, max.dist, by=signif(res(template)[1], 2))
}
if (!is.null(output_directory)) {
saveRDS(radii, file=paste0(output_directory, "\\PCF_radius.rds"))
}
if (mc.cores>1) {
cl <- snow::makeCluster(mc.cores, type="SOCK")
doSNOW::registerDoSNOW(cl)
pcf_brk <- foreach(i=yr.range, .combine=rbind, .inorder=FALSE, .packages=c("raster", "magrittr", "dplyr", "spatstat", "stringr", "maptools"), .export="create_raster") %dopar% {
add.df <- data.frame()
for (j in 1:12) {
brk_ind <- brks_only %>%
filter(year==i, month==j)
if (nrow(brk_ind)==0) next
brk_pts <- create_raster("no_cell", "brk", brk_ind, template)
if (!is.null(resample_dist)) {
brk_pts <- resample(brk_pts, resample_ras, "ngb")
}
brk_pts <- brk_pts %>%
rasterToPoints(spatial=TRUE) %>%
as.ppp.SpatialPointsDataFrame()
Window(brk_pts) <- window
if (is.na(level)) {
pcf <- pcf(brk_pts, r=radii, correction="translate")
g <- pcf$trans
} else {
env <- envelope(brk_pts, fun=pcf, nsim=simulations, r=radii, correction="translate")
g <- env$obs
}
g[is.infinite(g)] <- NA
ind.localmax <- which(diff(sign(diff(g)))==-2)+1
ind.firstpeak <- min(ind.localmax); ind.firstpeak <- ifelse(is.infinite(ind.firstpeak), NA, ind.firstpeak)
g.peak <- max(g[ind.localmax]); g.peak <- ifelse(is.infinite(g.peak), NA, g.peak)
ind.peak <- ifelse(is.na(g.peak), NA, which(g==g.peak))
ind.extrema <- c(which.min(g), which.max(g), ind.peak, ind.firstpeak)
add.df <- rbind(add.df, data.frame(year=i, month=j, type=c("min", "max", "peak", "first peak"), g=g[ind.extrema], dist=radii[ind.extrema],
upper.bound=ifelse(is.na(level), NA, env$hi[ind.extrema]), lower.bound=ifelse(is.na(level), NA, env$lo[ind.extrema]), nsim=simulations, sig.level=level))
if (!is.null(output_directory)) {
saveRDS(g, file=paste0(output_directory, "\\PCF_", i, "_", str_pad(j, width=2, side="left", pad="0"), ".rds"))
}
}
return(add.df)
}
snow::stopCluster(cl)
} else {
pcf_brk <- data.frame()
for (i in yr.range) {
for (j in 1:12) {
brk_ind <- brks_only %>%
filter(year==i, month==j)
if (nrow(brk_ind)==0) next
brk_pts <- create_raster("no_cell", "brk", brk_ind, template)
if (!is.null(resample_dist)) {
brk_pts <- resample(brk_pts, resample_ras, "ngb")
}
brk_pts <- brk_pts %>%
rasterToPoints(spatial=TRUE) %>%
as.ppp.SpatialPointsDataFrame()
Window(brk_pts) <- window
env <- envelope(brk_pts, fun=pcf, nsim=3, r=radii, correction="translate")
g <- env$obs; g[is.infinite(g)] <- NA
ind.localmax <- which(diff(sign(diff(g)))==-2)+1
ind.firstpeak <- min(ind.localmax); ind.firstpeak <- ifelse(is.infinite(ind.firstpeak), NA, ind.firstpeak)
g.peak <- max(g[ind.localmax]); g.peak <- ifelse(is.infinite(g.peak), NA, g.peak)
ind.peak <- ifelse(is.na(g.peak), NA, which(g==g.peak))
ind.extrema <- c(which.min(g), which.max(g), ind.peak, ind.firstpeak)
pcf_brk <- rbind(pcf_brk, data.frame(year=i, month=j, type=c("min", "max", "peak", "first peak"), g=g[ind.extrema], dist=radii[ind.extrema],
upper.bound=env$hi[ind.extrema], lower.bound=env$lo[ind.extrema], nsim=simulations, sig.level=level))
if (!is.null(output_directory)) {
saveRDS(g, file=paste0(output_directory, "\\PCF_", i, "_", str_pad(j, width=2, side="left", pad="0"), ".rds"))
}
}
}
}
return(list(df=pcf_brk, radius=radii))
}
jnghiem/bfasttools documentation built on Nov. 4, 2019, 3:02 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions pcf_calc
Documented in pcf_calc
#' Compute pair correlation function for breakpoint distribution
#'
#' This function computes the pair correlation function (PCF) for the point
#' pattern of breakpoints. The PCF estimates \emph{g} for a given radius
#' distance \emph{r}, where \emph{g} represents the degree of clustering.
#'
#' This function is a convenient wrapper for \code{spatstat::pcf()} and
#' \code{spatstat::envelope()}. Most default values in \code{pcf} and
#' \code{envelope} are accepted. However, the translation correction is used
#' here.
#'
#' If \code{level} is not set to \code{NA}, a simulation envelope will be
#' computed according to the number of simulations needed to achieve
#' \code{level}. However, the simulation process is very time-consuming so
#' setting \code{level=NA} may be best choice for reasonable processing times.
#'
#' The returned data frame has the following fields: \enumerate{ \item
#' \code{year}: The year of the observation. \item \code{month}: The month of
#' the observation. \item \code{type}: The type of statistic for \emph{g}. \item
#' \code{g}: The value of \emph{g} for the given \code{type}. \item \code{dist}:
#' The radius distance. \item \code{upper.bound, lower.bound}: The upper/lower
#' bound of the simulation envelope at distance \code{dist}. \item \code{nsim}:
#' The number of simulations to produce the envelope. \item \code{sig.level}:
#' The significance level of the pointwise Monte Carlo test. }
#'
#' @param bfast_in A data frame generated by \code{bfastSpatial()}.
#' @param level A number indicating the significance level for the pointwise
#' Monte Carlo test. See the documentation for \code{spatstat::envelope()} for
#' more details. If \code{level=NA}, no envelope will be simulated.
#' @param boundary A SpatialPolygons object defining the boundary of the PCF
#' analysis. If this argument is not supplied, all records in \code{bfast_in}
#' will be used.
#' @param template A Raster object with the same resolution, extent, and
#' projection as the raster dataset used to generate \code{bfast_in}.
#' Alternatively, an XML file generated by \code{create_raster_metadata()}.
#' @param resample_dist An optional number indicating the square dimension to
#' which the breakpoint pattern raster should be resampled. This argument is
#' in the same units as the projection of \code{template}. If not supplied, no
#' resampling will occur and points will be generated based on the resolution
#' of \code{template}.
#' @param max.dist The maximum radius distance to which the PCF calculation
#' should be performed. If not supplied, the default is the diagonal distance
#' across \code{boundary}.
#' @param output_directory An optional character file path to a directory to
#' which the PCF calculations for each year and month should be written in
#' \code{.rds} format. The radius distances for which the PCF was computed
#' will also be written here.
#' @param mc.cores A numeric indicating the number of cores to be used in
#' parallel computing.
#' @return A list of length 2. The first item is a data frame summarizing the
#' PCF results. The second item is the vector of radius distances.
#' @examples
#' \dontrun{
#' pcf_calc(bf_df, NA, rgdal::readOGR(dsn="C:/Desktop/shapefiles", layer="Mojave"), template, 2000, "C:/Desktop/PCF_files", 6)
#' }
#' @importFrom spatstat owin Window pcf envelope
#' @importFrom foreach %dopar% foreach
#' @importFrom maptools as.ppp.SpatialPointsDataFrame
#' @export
pcf_calc <- function(bfast_in, level=0.05, boundary, template, resample_dist=NULL, max.dist=NULL, output_directory=NULL, mc.cores=1) {
template <- create_template_raster(template)
if (!is.null(boundary)) {
cells <- extract(template, boundary, cellnumbers=TRUE, df=TRUE)[,"cell"]
} else {
cells <- 1:ncell(template)
}
results_cleaned <- bfast_in %>%
clean_bfast() %>%
filter(no_cell %in% cells, !is.na(start_date))
ext <- extent(boundary)
window <- owin(xrange=ext[1:2], yrange=ext[3:4])
brks_only <- filter(results_cleaned, brk==1)
years <- results_cleaned[,"year"]; yr.range <- min(years):max(years)
simulations <- round((2/level)-1)
if (is.null(max.dist)) {
max.dist <- ceiling(sqrt(diff(ext[1:2])^2+diff(ext[3:4])^2)/2)
}
if (!is.null(resample_dist)) {
resample_ras <- raster(ext, resolution=resample_dist, crs=proj4string(template))
radii <- seq(0, max.dist, by=signif(res(resample_ras)[1], 2))
} else {
radii <- seq(0, max.dist, by=signif(res(template)[1], 2))
}
if (!is.null(output_directory)) {
saveRDS(radii, file=paste0(output_directory, "\\PCF_radius.rds"))
}
if (mc.cores>1) {
cl <- snow::makeCluster(mc.cores, type="SOCK")
doSNOW::registerDoSNOW(cl)
pcf_brk <- foreach(i=yr.range, .combine=rbind, .inorder=FALSE, .packages=c("raster", "magrittr", "dplyr", "spatstat", "stringr", "maptools"), .export="create_raster") %dopar% {
add.df <- data.frame()
for (j in 1:12) {
brk_ind <- brks_only %>%
filter(year==i, month==j)
if (nrow(brk_ind)==0) next
brk_pts <- create_raster("no_cell", "brk", brk_ind, template)
if (!is.null(resample_dist)) {
brk_pts <- resample(brk_pts, resample_ras, "ngb")
}
brk_pts <- brk_pts %>%
rasterToPoints(spatial=TRUE) %>%
as.ppp.SpatialPointsDataFrame()
Window(brk_pts) <- window
if (is.na(level)) {
pcf <- pcf(brk_pts, r=radii, correction="translate")
g <- pcf$trans
} else {
env <- envelope(brk_pts, fun=pcf, nsim=simulations, r=radii, correction="translate")
g <- env$obs
}
g[is.infinite(g)] <- NA
ind.localmax <- which(diff(sign(diff(g)))==-2)+1
ind.firstpeak <- min(ind.localmax); ind.firstpeak <- ifelse(is.infinite(ind.firstpeak), NA, ind.firstpeak)
g.peak <- max(g[ind.localmax]); g.peak <- ifelse(is.infinite(g.peak), NA, g.peak)
ind.peak <- ifelse(is.na(g.peak), NA, which(g==g.peak))
ind.extrema <- c(which.min(g), which.max(g), ind.peak, ind.firstpeak)
add.df <- rbind(add.df, data.frame(year=i, month=j, type=c("min", "max", "peak", "first peak"), g=g[ind.extrema], dist=radii[ind.extrema],
upper.bound=ifelse(is.na(level), NA, env$hi[ind.extrema]), lower.bound=ifelse(is.na(level), NA, env$lo[ind.extrema]), nsim=simulations, sig.level=level))
if (!is.null(output_directory)) {
saveRDS(g, file=paste0(output_directory, "\\PCF_", i, "_", str_pad(j, width=2, side="left", pad="0"), ".rds"))
}
}
return(add.df)
}
snow::stopCluster(cl)
} else {
pcf_brk <- data.frame()
for (i in yr.range) {
for (j in 1:12) {
brk_ind <- brks_only %>%
filter(year==i, month==j)
if (nrow(brk_ind)==0) next
brk_pts <- create_raster("no_cell", "brk", brk_ind, template)
if (!is.null(resample_dist)) {
brk_pts <- resample(brk_pts, resample_ras, "ngb")
}
brk_pts <- brk_pts %>%
rasterToPoints(spatial=TRUE) %>%
as.ppp.SpatialPointsDataFrame()
Window(brk_pts) <- window
env <- envelope(brk_pts, fun=pcf, nsim=3, r=radii, correction="translate")
g <- env$obs; g[is.infinite(g)] <- NA
ind.localmax <- which(diff(sign(diff(g)))==-2)+1
ind.firstpeak <- min(ind.localmax); ind.firstpeak <- ifelse(is.infinite(ind.firstpeak), NA, ind.firstpeak)
g.peak <- max(g[ind.localmax]); g.peak <- ifelse(is.infinite(g.peak), NA, g.peak)
ind.peak <- ifelse(is.na(g.peak), NA, which(g==g.peak))
ind.extrema <- c(which.min(g), which.max(g), ind.peak, ind.firstpeak)
pcf_brk <- rbind(pcf_brk, data.frame(year=i, month=j, type=c("min", "max", "peak", "first peak"), g=g[ind.extrema], dist=radii[ind.extrema],
upper.bound=env$hi[ind.extrema], lower.bound=env$lo[ind.extrema], nsim=simulations, sig.level=level))
if (!is.null(output_directory)) {
saveRDS(g, file=paste0(output_directory, "\\PCF_", i, "_", str_pad(j, width=2, side="left", pad="0"), ".rds"))
}
}
}
}
return(list(df=pcf_brk, radius=radii))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.