R/lda_raster.R

Defines functions lda_raster

Documented in lda_raster

#' 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)

}
juoe/spatialtools documentation built on May 25, 2019, 6:25 p.m.