Commit 7b507821 authored by Monnet Jean-Matthieu's avatar Monnet Jean-Matthieu

Added error sampling in ABApredict, in log transform case

parent bf8dcb6e
......@@ -372,10 +372,11 @@ ABAmodelPlot <- function(modell, title=NULL, unit=NULL, disp.text=F, ...)
#' @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}
#' @param addError boolean. indicates whether errors sampled from a normal distribution N(0, sigma(residuals)) should be added to fitted values; implemented only for log transformation case
#' @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)
ABApredict <- function(modell, map.metrics, strata=NULL, addError=FALSE)
{
# create factor of strata if not existing
if (is.null(strata))
......@@ -401,6 +402,7 @@ ABApredict <- function(modell, map.metrics, strata=NULL)
dummy[[stratum]] <- raster::predict(map.metrics[[variables]], modell$model[[stratum]])
# back-transform
dummy[[stratum]] <- lidaRtRee::iBoxcoxTrBiasCor(dummy[[stratum]], modell$stats[stratum,"lambda"], modell$stats[stratum,"var.res"])
if (addError==TRUE) {warning("Error sampling not implemented in the boxcox transformation case")}
}
if (modell$stats[stratum,"transform"] == "log") # case of case of log-log transform
{
......@@ -409,13 +411,21 @@ ABApredict <- function(modell, map.metrics, strata=NULL)
newdata <- log(map.metrics[[variables]])
names(newdata) <- variables
dummy[[stratum]] <- raster::predict(newdata, modell$model[[stratum]])
dummy[[stratum]]<- exp(dummy[[stratum]]) * exp(modell$stats[stratum,"var.res"]/2)
if (addError==TRUE)
{
rastResidual <- dummy[[stratum]]
raster::values(rastResidual) <- stats::rnorm(length(rastResidual), 0, sqrt(modell$stats[stratum,"var.res"]))
dummy[[stratum]]<- exp(dummy[[stratum]] + rastResidual)
} else {
dummy[[stratum]]<- exp(dummy[[stratum]]) * exp(modell$stats[stratum,"var.res"]/2)
}
}
if (modell$stats[stratum,"transform"] == "none") # case of case of no transform
{
variables <- names(modell$model[[stratum]]$coefficients)
variables <- variables[variables !="(Intercept)"]
dummy[[stratum]] <- raster::predict(map.metrics[[variables]], modell$model[[stratum]])
if (addError==TRUE) {warning("Error sampling not implemented in the none transformation case")}
}
# set prediction outside of strata to NA
dummy[[stratum]][map.metrics$strata != stratum.ID] <- NA
......@@ -443,8 +453,8 @@ cleanPredictionRaster <- function(rast, minmax=c(-Inf, +Inf), mask=NULL)
rast <- rast * mask
}
#
rast[rast<minmax[1]] <- minmax[1]
rast[rast>minmax[2]] <- minmax[2]
rast[rast<=minmax[1]] <- minmax[1]
rast[rast>=minmax[2]] <- minmax[2]
#
return(rast)
}
......
......@@ -258,7 +258,9 @@ points2terrainStats <- function(p, centre=NULL, r=NULL)
#' dummy <- x$h[which(x$h>20 & x$h<30)]
#' data.frame(Tree.between.20.30=length(dummy)/area.ha, Tree.meanH=mean(dummy))
#' }
#' cloudTreeMetrics(llas, cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)), 10, res=0.5, func=user.func)
#' cloudTreeMetrics(llas,
#' cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)),
#' 10, res=0.5, func=user.func)
#' @export
#'
cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, func=stdTreeMetrics, ...)
......
......@@ -4,7 +4,7 @@
\alias{ABApredict}
\title{Mapping of ABA prediction models}
\usage{
ABApredict(modell, map.metrics, strata = NULL)
ABApredict(modell, map.metrics, strata = NULL, addError = FALSE)
}
\arguments{
\item{modell}{model returned by \code{\link{ABAmodel}} or \code{\link{ABAmodelCombineStrata}}}
......@@ -12,6 +12,8 @@ ABApredict(modell, map.metrics, strata = NULL)
\item{map.metrics}{RasterStack. metrics returned e.g by \code{\link[lidR]{grid_metrics}}}
\item{strata}{string. indicates which layer of map.metrics contains the \code{strata} in \code{modell}}
\item{addError}{boolean. indicates whether errors sampled from a normal distribution N(0, sigma(residuals)) should be added to fitted values; implemented only for log transformation case}
}
\value{
a raster of predictions obtained by applying the model \code{modell} to the observations in \code{map.metrics}
......
......@@ -52,7 +52,9 @@ user.func <- function(x, area.ha=NA)
dummy <- x$h[which(x$h>20 & x$h<30)]
data.frame(Tree.between.20.30=length(dummy)/area.ha, Tree.meanH=mean(dummy))
}
cloudTreeMetrics(llas, cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)), 10, res=0.5, func=user.func)
cloudTreeMetrics(llas,
cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)),
10, res=0.5, func=user.func)
}
\seealso{
\code{\link{treeSegmentation}}, \code{\link{treeExtraction}}, \code{\link{stdTreeMetrics}}
......
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