#' Mangrove Indices for Point-Centered Quarter Method
#'
#' This function allows you to calculate the Holdridge Complexity Index and Mean Stand Diameter for pcqm data.
#' @param samplingpoint Column name in data frame for sampling points (numerical). Default name is "samplingpoint". First sampling point must be 1.
#' @param dist Column name in data frame for distance from tree to sampling point. Default name is "dist". Values must be in meters.
#' @param species Column name in data frame for species. Default name is "species".
#' @param dbh Column name in data frame for diameter at breast height. Default name is "dbh". Values must be in centimeters.
#' @param height Column name in data frame for height (Optional). Default name is "height". Values must be in meters. When included, additional height-related outputs are displayed.
#' @param sizebin Logical argument for displaying size class plots (Default is FALSE). If sizebin=TRUE, plots of size class proportions are displayed for dbh (<5 cm dbh, 5–10 cm dbh, and >10 cm dbh) and height (<5 m, 5–10 m, and >10 m).
#' @keywords mangrove structure, pcqm, Holdridge Complexity Index, Mean Stand Diameter
#' @examples
#' pcqm.indices(mangrove_data)
#' pcqm.indices(mangrove_data, dbh = "Diameter", height = "Tree_Height", samplingpoint = "Sampling_Point", sizebin = T)
#' @export
pcqm.indices <- function(x,
samplingpoint = 'samplingpoint',
dist = 'dist',
species = 'species',
dbh = 'dbh',
height = 'height',
sizebin = F){
#Confirm it is a data frame
x <- as.data.frame(x)
#Load columns
x$SamplingPoint <- x[,samplingpoint]
x$Distance <- x[,dist]
x$Species <- x[,species]
x$dbh <- x[,dbh]
if("height" %in% colnames(x)) x$height <- x[,height]
## Get the functions from the plyr namespace
summarize = get("summarize", asNamespace('plyr'))
transform = get("transform", asNamespace('plyr'))
#Check for NA cells
if(any(is.na(x)) == TRUE) stop("Data frame cannot contain missing values (NAs).")
# Pre calculations for later use
spp <- length(unique(x$Species))
x$count <- 1:nrow(x)
cat("\n ---------------------------------------------------")
cat("\n Point-Centered Quarter Method Indices")
cat("\n ---------------------------------------------------\n\n")
# Calculate Density
x$static <- 1
densum <- plyr::ddply(x, "static", transform, sum.n = max(count), meandist= mean(Distance), stemcalc=((1/(mean(Distance)^2))*1000))
densum1 <- plyr::ddply(densum, "Species", summarize, Number = length(static), Proportion = (length(static) / max(sum.n)), stemcalc=max(stemcalc))
densum2 <- plyr::ddply(densum1, "Species", transform, Stems_Per_0.1_Ha = round(Proportion * stemcalc))
totden <- sum(densum2$Stems_Per_0.1_Ha)
# Calculate Basal Area
densum$BA <- (3.14159265359 * (densum$dbh/2) ^ 2) / 10000
densumBA <- plyr::ddply(densum, "Species", summarize, MeanBA = (mean(BA)))
densumBA <- plyr::join(densum2, densumBA, by="Species")
densumBA$Basal_Area <- densumBA$Stems_Per_0.1_Ha * densumBA$MeanBA
totbas <- sum(densumBA$Basal_Area)
# Get height of three tallest trees
if("height" %in% colnames(x)) x <- x[order(-x$height),]
if("height" %in% colnames(x)) x2 <- head(x,3)
if("height" %in% colnames(x)) heightcalc <- mean(x2$height)
# Calculate Holdridge Complexity Index
if("height" %in% colnames(x)) CI <- round((heightcalc * totbas * totden * spp) / 1000, digits=2)
# Print the Output from HCI
if("height" %in% colnames(x)) cat(" Complexity Index\n -----\n")
if("height" %in% colnames(x)) cat(paste(" ", CI))
if("height" %in% colnames(x)) cat("\n\n")
# Calculate Mean Stand Diameter
MSD <- round(sqrt((totbas * 12732.39)/totden), digits=2)
# Print output from MSD
cat(" Mean Stand Diameter\n -----\n")
cat(paste(" ", MSD,"\n\n"))
if(sizebin==T){
# Plot of dbh proportions
x$size <- "<5 cm"
x$size[x$dbh >=5 & x$dbh<=10] <- "5–10 cm"
x$size[x$dbh > 10] <- ">10 cm"
x$sizefac <- factor(x$size, levels= c("<5 cm", "5–10 cm", ">10 cm"))
plotsize <- table(x$sizefac)
par(ask = TRUE)
barplot(plotsize / sum(plotsize), ylab="Proportion of forest trees", xlab="Tree diameter size class", cex.lab = 1.2, cex.names=1.1, cex.axis = 1.1, ylim=c(0,1), col = c("burlywood3", "tan3", "saddlebrown"))
# Plot of height proportions
if("height" %in% colnames(x)) x$size2 <- "<5 m"
if("height" %in% colnames(x)) x$size2[x$height >=5 & x$height<=10] <- "5–10 m"
if("height" %in% colnames(x)) x$size2[x$height > 10] <- ">10 m"
if("height" %in% colnames(x)) x$sizefac2 <- factor(x$size2, levels= c("<5 m", "5–10 m", ">10 m"))
if("height" %in% colnames(x)) plotsize2 <- table(x$sizefac2)
if("height" %in% colnames(x)) barplot(plotsize2 / sum(plotsize2), ylab="Proportion of forest trees", xlab="Tree height size class", cex.lab = 1.2, cex.names=1.1, cex.axis = 1.1, ylim=c(0,1), col = c("darkolivegreen3", "forestgreen", "darkgreen"))
#Scatterplot of dbh vs height
if("height" %in% colnames(x)) plot(x$dbh, x$height, pch=21, xlab="Tree DBH (cm)", ylab="Tree Height (m)", cex.lab = 1.2, cex.axis = 1.1, col="black", bg = c("snow", "gray75", "black")[x$sizefac], ylim=c(0,max(x$height)+1), xlim=c(0,max(x$dbh)+1), yaxs="i", xaxs="i")
if("height" %in% colnames(x)) legend("bottomright",legend=c("<5 cm dbh","5–10 cm dbh",">10 cm dbh"),pch=c(21,21,21),col = "black", pt.bg = c("snow", "gray75", "black"), cex=1.2, pt.cex=1.2, title="DBH Size Class")
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.