#' Apply LDA to raster stack
#'
#' @param x Raster stack of predictors
#' @param y Raster layer with dependent categorical variable
#' @param nsamples Number of samples to use for calculation of LDA
#' @param ncores Number of cores to use
#' @param filename Optional. filename for writing raster to disk.
#' @param ... arguments passed on to PCA function
#'
#' @return A list with 1) the raster stack with the LDA components and 2) the LDA object
#' @export
lda_raster <- function(x, y, nsamples = 200, ncores, filename = NULL, ...) {
maparea <- as.numeric(table(raster::values(y)[which(!is.na(raster::values(y)))]))
priors <- maparea/sum(maparea)
nlds <- length(priors)-1
s <- raster::sampleStratified(y, nsamples, sp = TRUE)
ext <- extract_parallel(x, s, ncores)
ext <- ext[complete.cases(ext),]
colnames(ext) <- names(x)
ext <- scale(ext)
dat <- cbind(data.frame(y = as.factor(s[[2]])), ext)
l <- MASS::lda(y~., dat, prior = priors)
lda_stack <- raster::stack(x)
bs <- raster::blockSize(lda_stack)
cl <- parallel::makeCluster(ncores)
doParallel::registerDoParallel(cl)
values <- foreach::foreach(i=1:bs$n, .packages=c("raster", "MASS"), .combine='rbind') %dopar% {
v <- getValues(x, row=bs$row[i], nrows=bs$nrows[i])
m <- matrix(NA, nrow(v), nlds)
newdata <- as.data.frame(v[which(rowSums(is.na(v))==0),])
preds <- predict(l, newdata)[["x"]]
m[which(rowSums(is.na(v))==0),] <- preds
m
}
parallel::stopCluster(cl)
lda_stack <- setValues(lda_stack, values)
names(lda_stack) <- paste0("LD", 1:nlayers(lda_stack))
if(!is.null(filename)) {
writeRaster(lda_stack, filename, overwrite=TRUE)
}
list(lda_stack, l)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.