R/pca_raster.R

In juoe/spatialtoolbox: Collection of functions for doing spatial stuff.

#' Apply PCA to raster stack
#'
#' @param x Raster stack
#' @param nsamples Number of samples to use for calculation of PCA
#' @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 PCA components and 2) the PCA object
#' @export

pca_raster <- function(x, nsamples = 10000, ncores, filename = NULL, ...) {

  s <- xyFromCell(x, sample(1:ncell(x), nsamples), spatial = TRUE)

  ext <- extract_parallel(x, s, ncores)
  ext <- ext[complete.cases(ext),]

  p <- prcomp(ext, ...)

  pca_stack <- predict_raster(x, p, ncores, stack = TRUE)
  names(pca_stack) <- paste0("PC", 1:nlayers(pca_stack))

  if(!is.null(filename)) {
    writeRaster(pca_stack, filename, overwrite=TRUE)
  }

  list(pca_stack, p)

}