Commit 9a500ebd authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu

Update ABAmodel to comply with raster outputs of grid_metrics

parent 00de3c46
......@@ -305,11 +305,11 @@ ABAmodelCombineStrata <- function(modells, plotsId=NULL)
return(modell)
}
####################################
#' plots observed and predicted values of a model returned by \code{\link{ABA.model}}
#' plots observed and predicted values of a model returned by \code{\link{ABAmodel}}
#'
#' ADD DESCRIPTION
#'
#' @param modell list .as returned by \code{\link{ABA.model}}
#' @param modell list. as returned by \code{\link{ABAmodel}}
#' @param title string. the plot title
#' @param unit numeric. the plot title (between parenthesis)
#' @param disp.text boolean indicates if points should be labelled with id
......@@ -346,51 +346,53 @@ ABAmodelPlot <- function(modell, title=NULL, unit=NULL, disp.text=F)
####################################
#' Mapping of ABA prediction models
#'
#' applies calibrated area-based prediction models output of \code{\link{ABA.model}} to a dataframe of new metrics
#' applies calibrated area-based prediction models output of \code{\link{ABAmodel}} to a RasterStack of new metrics
#'
#' @param modell a model returned by \code{\link{ABAmodel}} or \code{\link{ABAmodelCombineStrata}}
#' @param map.metrics metrics returned e.g by \code{\link[lidR]{grid_metrics}}
#' @param strata a factor vector. levels corresponding to \code{strata} in \code{modell}
#' @return a vector of predictions obtained by applying the model \code{modell} to the observations (lines) in \code{map.metrics}
#' @param modell model returned by \code{\link{ABAmodel}} or \code{\link{ABAmodelCombineStrata}}
#' @param map.metrics RasterStack. metrics returned e.g by \code{\link[lidR]{grid_metrics}}
#' @param strata string. indicates which layer of map.metrics contains the \code{strata} in \code{modell}
#' @return a raster of predictions obtained by applying the model \code{modell} to the observations in \code{map.metrics}
#' @export
#'
ABApredict <- function(modell, map.metrics, strata=NULL)
{
map.metrics <- as.data.frame(map.metrics)
# create factor of strata
# create factor of strata if not existing
if (is.null(strata))
{
strata <- as.factor(rep("1",nrow(map.metrics)))
row.names(modell$stats) <- "1"
map.metrics$strata <- 1
map.metrics$strata <- as.factor(map.metrics$strata)
levels(map.metrics$strata) <- data.frame(ID=1, propriete="1")
}
# create final vector
dummy <- rep(1.1, nrow(map.metrics))
dummy[1:nrow(map.metrics)] <- NA
#
dummy <- list()
# loop on strata
for (i in levels(strata))
for (i in 1:nrow(levels(map.metrics[["strata"]])[[1]]))
{
# extract stratum
stratum <- which(strata==i)
# predict
if (modell$stats[i,"transform"] == "boxcox") # case of Box-Cox transform
stratum <- as.character(levels(map.metrics[["strata"]])[[1]][i, 2])
stratum.ID <- levels(map.metrics[["strata"]])[[1]][i, "ID"]
# predict on all cells
if (modell$stats[stratum,"transform"] == "boxcox") # case of Box-Cox transform
{
# apply linear model
dummy[stratum] <- stats::predict.lm(modell$model[[i]], newdata=map.metrics[stratum,])
dummy[[stratum]] <- raster::predict(map.metrics[[names(modell$model[[stratum]]$coefficients)[-1]]], modell$model[[stratum]])
# back-transform
dummy[stratum] <- iBoxcoxTrBiasCor(dummy[stratum], modell$stats[i,"lambda"], modell$stats[i,"var.res"])
dummy[[stratum]] <- lidaRtRee::iBoxcoxTrBiasCor(dummy[[stratum]], modell$stats[stratum,"lambda"], modell$stats[i,"var.res"])
}
if (modell$stats[i,"transform"] == "log") # case of case of log-log transform
if (modell$stats[stratum,"transform"] == "log") # case of case of log-log transform
{
dummy[stratum] <- stats::predict.lm(modell$model[[i]], newdata=log(map.metrics[stratum,]))
dummy[stratum] <- exp(dummy[stratum]) * exp(modell$stats[i,"var.res"]/2)
dummy[[stratum]] <- raster::predict(log(map.metrics[[names(modell$model[[stratum]]$coefficients)[-1]]]), modell$model[[stratum]])
dummy[[stratum]]<- exp(dummy[[stratum]]) * exp(modell$stats[stratum,"var.res"]/2)
}
if (modell$stats[i,"transform"] == "none") # case of case of no transform
if (modell$stats[stratum,"transform"] == "none") # case of case of no transform
{
dummy[stratum] <- stats::predict.lm(modell$model[[i]], newdata=map.metrics[stratum,])
dummy[[stratum]] <- raster::predict(map.metrics[[names(modell$model[[stratum]]$coefficients)[-1]]], modell$model[[stratum]])
}
# set prediction outside of strata to NA
dummy[[stratum]][map.metrics$strata != stratum.ID] <- NA
}
return(dummy)
# merge strata results
return(Reduce(function(...)raster::merge(...,tolerance=1),dummy))
}
####################################
#' applies thresholds and mask to a prediction map
......@@ -419,7 +421,7 @@ cleanPredictionRaster <- function(r, minmax=c(-Inf, +Inf), mask=NULL)
}
####################################
#' computes inference from area-based model and predicted values
#' @param modell a model returned by by \code{\link{ABA.model}} or \code{\link{combine.stratified.ABA.models}}
#' @param modell a model returned by by \code{\link{ABAmodel}} or\code{\link{ABAmodelCombineStrata}}
#' @param r.predictions raster of predicted values
#' @param type string vector specifying which estimators should be computed (one or several in "SRS", "ED", "D", "STR", "SYNT")
#' @param r.mask raster to mask region of interest (NA values), may contain post-stratification categories (should be integer, positive values)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment