#' get_patches
#'
#' Find hot and cold patches in a numeric matrix or raster, and calculate patch statistics.
#' @param img A numeric temperature matrix (such as that returned from
#' \code{Thermimage::}\code{\link[Thermimage]{raw2temp}}) or raster.
#' @param id The image ID (optional). Useful when iterating over
#' numerous images.
#' @param style Style to use when calculating neighbourhood weights using
#' \code{spdep::}\code{\link[spdep]{nb2listw}}. Defaults to 'C' (globally
#' standardised).
#' @param img_proj Spatial projection. Optional, but necessary for geographic
#' data to plot correctly.
#' @param img_extent Spatial extent. Optional, but necessary for geographic
#' data to plot correctly.
#' @param return_vals Which values to return? Any combination of the dataframe
#' (\code{df}), SpatialPolygonsDataFrame of hot and cold patches
#' (\code{patches}) and patch statistics dataframe (\code{pstats}).
#' @return A list containing:
#' \item{df}{A dataframe with one row for each pixel, and variables denoting:
#' the pixel value (val); the spatial location of the pixel (x and y);
#' its patch classification (G_bin) into a hot (1), cold (-1) or no patch (0)
#' according to the Z value (see \code{spdep::}\code{\link[spdep]{localG}});
#' the unique ID of the patch in which the pixel fell;
#' the image ID (if applicable); and the original spatial location of the pixel
#' (x_orig and y_orig).}
#' \item{patches}{A SpatialPolygonsDataFrame of hot and cold patches, named
#' according to \code{id} (if applicable). Hot patches have a value of 1,
#' and cold patches a value of -1.}
#' \item{pstats}{A dataframe with patch statistics for hot patches and cold
#' patches, respectively. See \code{\link{patch_stats}} for details of all the
#' statistics returned.}
#' @importFrom rlang .data
#' @examples
#'
#' \dontrun{
#' # FLIR temperature matrix ----------------------------------------
#' # Find hot and cold patches
#' flir_results <-
#' get_patches(img = flir11835$flir_matrix,
#' id = flir11835$photo_no)
#'
#' # Look at the results for individual pixels
#' head(flir_results$df)
#'
#' # Look at the patch statistics for hot and cold patches
#' flir_results$pstats
#'
#' # Plot the patches
#' sp::plot(flir_results$patches)
#'
#' # Plot using ThermStats::plot_patches
#' plot_patches(df = flir_results$df,
#' patches = flir_results$patches,
#' print_plot = TRUE,
#' save_plot = FALSE)
#'
#' # Plot using landscapemetrics package
#' flir_patch_rast <- raster::rasterFromXYZ(flir_results$df[,c("x", "y", "G_bin")])
#' if (requireNamespace("landscapemetrics", quietly = TRUE)){
#' show_landscape(flir_patch_rast)
#' show_patches(flir_patch_rast, class = "all", labels = FALSE)
#' show_cores(flir_patch_rast, class = c(1, -1), labels = FALSE)
#' }
#'
#' # Worldclim2 temperature raster ----------------------------------
#' # Dataset 'sulawesi_temp' represents mean January temperature for
#' # the island of Sulawesi
#'
#' # Define projection and extent
#' img_proj <- raster::projection(sulawesi_temp)
#' img_extent <- raster::extent(sulawesi_temp)
#'
#' # Find hot and cold patches
#' worldclim_results <-
#' get_patches(img = sulawesi_temp,
#' id = "sulawesi",
#' img_proj = img_proj,
#' img_extent = img_extent)
#'
#' # Look at the results for individual pixels
#' head(worldclim_results$df)
#'
#' # Look at the patch statistics for hot and cold patches
#' worldclim_results$pstats
#'
#' # Plot the patches
#' sp::plot(worldclim_results$patches)
#'
#' # Plot using ThermStats::plot_patches
#' plot_patches(df = worldclim_results$df,
#' patches = worldclim_results$patches,
#' print_plot = TRUE,
#' save_plot = FALSE)
#' }
#'
#' @importClassesFrom sp SpatialPolygonsDataFrame SpatialPointsDataFrame
#' @export
#'
get_patches <- function(img,
id = NULL,
style = "C",
img_proj = NULL,
img_extent = NULL,
return_vals = c("df","patches","pstats")) {
# Setup --------------------------------------------------------------------
message("Getting patches")
# Get dimensions
nrows <- nrow(img)
ncols <- ncol(img)
# Determine whether matrix or raster
if(is.matrix(img)){
# Melt to dataframe
df <- reshape2::melt(img,
varnames = c("y", "x"),
value.name = "val")
}else if(class(img)[1] == "RasterLayer"){
# Coerce to dataframe
df <- raster::as.data.frame(img, xy = TRUE)
colnames(df)[3] <- "val"
}
# Remove NA values
df <- df[!(is.na(df[, "val"])),]
# Neighbour weights -------------------------------------------------------
message("\t...calculating neighbourhood weights")
# Identify nearest 8 neighbours
nr_neigh <- spdep::knearneigh(cbind(df$x, df$y), k = 8)
nr_neigh <- spdep::knn2nb(nr_neigh)
# Include self to calculate G* variant of local G statistic
nr_neigh <- spdep::include.self(nr_neigh)
# Get neighbour weights
nb_weights <- spdep::nb2listw(nr_neigh,
style = style,
zero.policy = FALSE)
# Local Getis-Ord ---------------------------------------------------------
message("\t...calculating local G* statistic")
local_g <-
spdep::localG(x = spdep::spNamedVec("val", df),
listw = nb_weights,
zero.policy = FALSE,
spChk = NULL)
# Add Z-values to dataframe
df <- data.frame(df, Z_val = matrix(local_g))
rm(nr_neigh, nb_weights, local_g)
# Define significance categories, using critical values defined in Getis and
# Ord (1996 )
get_critical_Z <- function(n){
if(n >= 1 & n < 10){
critical_z <- 1.6450
}else if(n >= 10 & n < 30){
critical_z <- 2.5683
}else if(n >= 30 & n < 50){
critical_z <- 2.9291
}else if(n >= 50 & n < 100){
critical_z <- 3.0833
}else if(n >= 100 & n < 1000){
critical_z <- 3.2889
}else if(n > 1000){
critical_z <- 3.8855
}else{
critical_z <- NA
}
}
critical_Z <- get_critical_Z(nrow(df))
df$G_bin <- ifelse(df$Z_val >= critical_Z, 1,
ifelse(df$Z_val <= -critical_Z, -1, 0))
# Patch statistics --------------------------------------------------------
message("\t...matching pixels to hot and cold patches")
# 1. Create layer with one polygon for each hot/cold patch
# Dataframe to matrix
patch_mat <- reshape2::acast(df, y ~ x, value.var = "G_bin")
# Flip
patch_mat <-
Thermimage::mirror.matrix(Thermimage::rotate180.matrix(patch_mat))
# Rasterise matrix
# patches <- raster::raster(patch_mat)
patches <-
raster::raster(patch_mat,
xmn=0, xmx=ncols,
ymn=0, ymx=nrows)
# Specify coordinates and img_projection (if applicable)
if(!(is.null(img_proj))){
raster::extent(patches) <- img_extent
raster::projection(patches) <- img_proj
}
# Raster to dissolved polygons
patches <- suppressMessages(raster::rasterToPolygons(patches, dissolve = TRUE))
patches <- raster::disaggregate(patches)
# 2. Assign to each temperature cell the ID of the patch that it falls into
# Matrix to raster to points
img <- reshape2::acast(df, y ~ x, value.var = "val")
# Flip
img <-
Thermimage::mirror.matrix(Thermimage::rotate180.matrix(img))
# Matrix to raster
raw <- raster::raster(img,
xmn=0, xmx=ncols,
ymn=0, ymx=nrows)
# Specify coordinates and img_projection (if applicable)
if(!(is.null(img_proj))){
raster::extent(raw) <- img_extent
raster::projection(raw) <- img_proj
}
# Raster to points
raw <- raster::rasterToPoints(raw, spatial = TRUE)
# Return overlay list where each element is the point, and within that the
# value is the hot/coldspot classification and the row name is the patch ID
n_features <- length(raw)
# Define function to overlay in subsets if too many features
over_sub <- function(min_i, max_i){
subdat <- raw[(min_i + 1):max_i,]
patchID <- sp::over(subdat, patches, returnList = TRUE)
return(patchID)
}
# Define function to reformat each dataframe within the list
reform <- function(x) {
x$patchID <- row.names(x)
colnames(x) <- c("G_bin", "patchID")
return(x)
}
if(n_features > 10^5){
lims <- seq(0, n_features, 10^5)
lims[lims == utils::tail(lims, n = 1)] <- n_features
# Overlay by subset
patchID <-
lapply(1:(length(lims) - 1), function(i){
over_sub(lims[i], lims[i + 1])
})
patchID <-
lapply(patchID, function(x){
# Reform each element (pixel) into a dataframe
df <-
lapply(x, function(y){
reform(y)
})
# Bind pixel dataframes within subsets
df <- do.call("rbind", df)
return(df)
})
# Bind all subsets
patchID <- do.call("rbind", patchID)
}else{
patchID <- sp::over(raw, patches, returnList = TRUE)
# Reform each element (pixel) into a dataframe
patchID <- lapply(patchID, reform)
# Bind pixel dataframes
patchID <- do.call("rbind", patchID)
}
# Recreate dataframe
df <- data.frame(raw)[-4]
colnames(df)[which(names(df) == "layer")] <- "val"
df[, c("G_bin", "patchID")] <- patchID
rm(raw, patchID)
# Function to determine what the first character of a string is
char_is_numeric <- function(str){
str <- as.character(str)
firstchar <- unlist(strsplit(str, ""))[1]
firstchar <- suppressWarnings(as.numeric(firstchar))
is_numeric <- ifelse(is.na(firstchar), FALSE, TRUE)
return(is_numeric)
}
if(!(is.null(id))){
if(char_is_numeric(id)) id <- as.character(paste("X", id, sep = ""))
df$id <- id
names(patches) <- as.character(id)
}
### 3. Calculate patch stats
message("\t...calculating patch statistics")
# Calculate median for each patch --------------------------------------------
if (requireNamespace("dplyr", quietly = TRUE)) {
patch_val <-
dplyr::summarise(dplyr::group_by(df[df$G_bin != 0, ],
patchID, .data$G_bin),
val = stats::median(.data$val, na.rm = TRUE))
patch_val <- as.data.frame(patch_val)
}else{
patchIDs <- unique(df[df[,"G_bin"] != 0, "patchID"])
patch_val <-
lapply(patchIDs, function(x){
patchID <- x
G_bin <- df[df[,"patchID"] == x,"G_bin"][1]
val <- stats::median(df[df[,"patchID"] == x,"val"], na.rm = TRUE)
return(data.frame(patchID = patchID,
G_bin = G_bin,
val = val))
})
patch_val <- do.call("rbind", patch_val)
}
# Calculate stats for each class --------------------------------------------
pstats <-
lapply(c(1,-1), function(class){
# Spatial stats
results <- patch_stats(mat = patch_mat, class = class)
# Patch number
results[, "abundance"] <-
length(unique(df[df[,"G_bin"] == class, "patchID"]))
# Patch density
results[, "density"] <-
results[, "abundance"] / results[, "total_area"]
# Summary stats across patches (based on median values)
if(length(patch_val[patch_val[,"G_bin"] == class, "val"]) > 0){
sumstats <-
t(as.matrix(summary(patch_val[patch_val[,"G_bin"] == class, "val"])))
cols <- tolower(gsub(" ", "", gsub("[.]", "", colnames(sumstats))))
cols[grep("1stqu", cols)] <- "lowerquant"
cols[grep("3rdqu", cols)] <- "upperquant"
colnames(sumstats) <- cols
# If no patches in this class, fill manually
}else{
sumstats <- data.frame(rbind(rep(NA,6)))
colnames(sumstats) <-
c("min","lowerquant","median","mean","upperquant","max")
}
results <- cbind(results, sumstats)
# Remove class column
results <- results[-1]
# Add patch class prefix to column names
prefix <- ifelse(class == 1, "hot", "cold")
colnames(results) <- paste(prefix, colnames(results), sep = "_")
# Return
return(results)
})
# Bind class results together
pstats <- do.call("cbind", pstats)
# Add image ID if available
if(!(is.null(id))) pstats[,"id"] <- id
# Return results ------------------------------------------------------------
if(length(return_vals) == 1){
return_vals <- eval(parse(text = return_vals))
}else{
return_vals <- sapply(return_vals, function(x) eval(parse(text = x)))
}
return(return_vals)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.