Nothing
if (getRversion() >= "2.15.1") {
utils::globalVariables(c(
"indice_line", "indice_col", "V1"
))
}
#' Summarize by plot (or subplot) the posterior distribution of AGB values
#'
#' @description
#' This function summarizes the matrix `AGB_val` given by the function [AGBmonteCarlo()] by plot. Or just do the sums
#' for each plot of the AGB if the argument `AGB_val` is the resulting vector from the function [computeAGB()].
#'
#' @details
#' If some trees belong to an unknown plot (i.e. NA value in the plot arguments), their AGB values are randomly assigned
#' to a plot at each iteration of the AGB monte Carlo approach. Or discarded when using output from [computeAGB()].
#'
#' The `drawPlot` argument is a logical that if it is set `TRUE`, a graph will appear with the plot given on absciss and the value
#' of AGB on ordinate, the red segments are the quantile, if `AGB_val` is the result of the function [AGBmonteCarlo()].
#' If the `subplot` arguments is set and the `drawPlot` is set `TRUE`, a graph is drawn with the spatialisation of the plots.
#'
#' @param AGB_val Matrix resulting from the function [AGBmonteCarlo()] (AGB_val element of the list),
#' or just the output of the function [AGBmonteCarlo()]. Or the output of the function [computeAGB()]
#' @param plot Vector with the code of plot
#' @param drawPlot a logical to draw the plot (see Details)
#' @param subplot Data frame, output of the function [cutPlot()]
#'
#' @return a data frame where:
#' - `plot`: the code of the plot
#' - `AGB`: AGB value at the plot level
#' - `Cred_2.5`: the quantile 2.5\% for the plot (when output of [AGBmonteCarlo()] is used)
#' - `Cred_97.5`: the quantile 97.5\% for the plot (when output of [AGBmonteCarlo()] is used)
#'
#' If the `subplot` is set, the output is a list with the previous data frame and a simple features (sf) geometry object.
#'
#'
#' @export
#'
#' @importFrom data.table data.table := first setDT
#' @importFrom grDevices terrain.colors
#' @importFrom graphics segments
#' @importFrom stats quantile
#' @examples
#'
#' # Load a database
#' data(NouraguesHD)
#' data(KarnatakaForest)
#'
#' # Modelling height-diameter relationship
#' HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, method = "log2")
#'
#' # Retrieving wood density values
#' \donttest{
#' KarnatakaWD <- getWoodDensity(KarnatakaForest$genus, KarnatakaForest$species,
#' stand = KarnatakaForest$plotId
#' )
#' }
#'
#' # Propagating errors
#' \donttest{
#' filt <- KarnatakaForest$plotId %in% c("BSP20", "BSP14")
#' resultMC <- AGBmonteCarlo(
#' D = KarnatakaForest$D[filt], WD = KarnatakaWD$meanWD[filt],
#' errWD = KarnatakaWD$sdWD[filt], HDmodel = HDmodel
#' )
#'
#' plot <- KarnatakaForest$plotId[ filt ]
#'
#' # The summary by plot
#' summaryByPlot(AGB_val = resultMC$AGB_simu, plot)
#'
#' # The summary by plot for computeAGB
#' H <- retrieveH(KarnatakaForest$D[filt], model = HDmodel)$H
#' AGB <- computeAGB(KarnatakaForest$D[filt], WD = KarnatakaWD$meanWD[filt], H = H)
#' summaryByPlot(AGB, plot)
#' }
summaryByPlot <- function(AGB_val, plot, drawPlot = FALSE, subplot = NULL) {
# parameters verification -------------------------------------------------
if (is.list(AGB_val)) {
AGB_val <- AGB_val$AGB_simu
}
if (!is.matrix(AGB_val) && !is.vector(AGB_val)) {
stop(
"The AGB_val must be a matrix you have for the result of the function ",
"'AGBmonteCarlo', or just the result of the function. ",
"Or the result from the function 'computeAGB'"
)
}
if (length(plot) != ifelse(is.matrix(AGB_val), nrow(AGB_val), length(AGB_val))) {
stop("Your 'plot' vector have not the same length as your number of row in the matrix")
}
if (!is.null(subplot)) {
if (!requireNamespace("sf", quietly = TRUE)) {
warning(
'To use this part of the function you must have the "sf" library\n\n',
'\t\tinstall.packages("sf")'
)
subplot <- NULL
}
if (!any(subplot$subplot %in% plot)) {
warning(
"The subplot parameter do not correspond to any plot"
)
subplot <- NULL
}
}
# function if it's a vector -----------------------------------------------
if (is.vector(AGB_val)) {
data <- data.table(AGB = AGB_val, plot = plot)
data <- na.omit(data)
AGB <- data[, .(AGB = sum(AGB)), by = plot]
} else {
# function if it's a matrix -----------------------------------------------
Plot <- data.table(plot = plot)
indice_tree <- Plot[is.na(plot), .I, by = plot]
# filter if there is there is NA in the AGB_val
filter <- rowSums(is.na(AGB_val)) > 0
# take the first tree in the database by plot
indice_first <- Plot[!is.na(plot) & !filter, .(indice_line = first(.I), plot = unique(plot)), by = plot]
# if there is trees without label -----------------------------------------
if (nrow(indice_tree) != 0) {
# Create a table with I : indice of the tree to distribute
# indice_line : a random sample of matrix line inside the plot without the trees to distribute
# indice_col : the column index of the matrix by I
mySample <- function(plot1, n) {
return(samples = indice_first[, sample(indice_line, n, replace = TRUE)])
}
n <- ncol(AGB_val)
samples <- indice_tree[, .(indice_line = mySample(plot, n), indice_col = 1:n), by = I]
# remove the index for the tree to distribute when it is NA
samples <- samples[ !(I %in% filter) ]
}
# function summary --------------------------------------------------------
mySummary <- function(x, matrix) {
# deal with plots with only one tree individual
if(!is.null(nrow(matrix[x, ]))){
resAGB <- colSums(matrix[x, ], na.rm = TRUE)
}else{
resAGB <-matrix[x, ]
}
# if there is trees without label
if (nrow(indice_tree) != 0) {
subsample <- samples[x[1] == indice_line, ]
# if the tree belong among the subplot
if (nrow(subsample) != 0) {
# sum for the trees I whose are not in the subplot, by column
sums <- subsample[, sum(matrix[I, unique(indice_col)], na.rm = TRUE), by = indice_col]
sums[is.na(V1), V1 := 0] # if there is any NA update the table
resAGB[sums$indice_col] <- resAGB[sums$indice_col] + sums$V1 # sum the result of the table
}
}
return(list(
AGB = mean(resAGB),
Cred_2.5 = quantile(resAGB, probs = 0.025),
Cred_97.5 = quantile(resAGB, probs = 0.975)
))
}
AGB <- Plot[!is.na(plot), mySummary(.I, AGB_val), by = plot]
if (drawPlot) {
with(AGB[order(AGB)], {
plot(AGB,
pch = 20, xlab = "", ylab = "AGB (Mg/ha)", ylim = range(Cred_2.5, Cred_97.5),
las = 1, cex.lab = 1.3, xaxt = "n", main = "AGB by plot"
)
axis(1, at = seq(length(AGB)), labels = plot, las = 2)
segments(x0 = seq(length(AGB)), y0 = Cred_2.5, y1 = Cred_97.5, col = "red")
})
}
}
if (!is.null(subplot)) {
setDT(subplot)
list_poly <- lapply(split(subplot, by = "subplot"), function(data) {
mat <- data[order(corner), .(XAbs, YAbs)]
mat <- as.matrix(rbind(mat, mat[1, ]))
output <- list()
output$poly <- sf::st_polygon(list(mat))
output$AGB <- AGB[ plot == unique(data$subplot), AGB] * data[, 10000 / (diff(range(XRel)) * diff(range(YRel)))]
output
})
sf_obj <- sf::st_sf(polygon = lapply(list_poly, "[[", 1), plot = names(list_poly), AGB = sapply(list_poly, "[[", 2))
if (drawPlot) {
plot(sf_obj["AGB"], main = "AGB (Mg / ha)", key.pos = 1, pal = function(n) {
rev(terrain.colors(n))
})
}
return(list(AGB = as.data.frame(AGB), polygon = sf_obj))
}
return(as.data.frame(AGB))
}
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.