R/phenexRaster.R

In joaofgoncalves/phenexRaster: Apply phenex analysis to raster objects

#' Get matrix with phenex modelled values
#'
#' Just an ancillary function to push data from \code{NDVI} objects generated by
#' \code{\link[phenex]{modelNDVI}} and construct a data matrix
#'
#' @param ndviObj An \code{NDVI} object generated by \code{\link[phenex]{modelNDVI}}
#'
#' @return A matrix containing modelled values

getModelledValuesToMatrix <- function(ndviObj){

  nr <- length(ndviObj)
  nc <- length(ndviObj[[1]]@modelledValues)

  ndviMat <- matrix(NA, nrow = nr, ncol = nc)

  for(i in 1:nr){
    ndviMat[i,] <- ndviObj[[i]]@modelledValues
  }

  return(ndviMat)

}

#' Applies phenex functions to a raster object
#'
#' A set of ancillary function used to apply phenex functionalities to a \code{Raster*} object.
#' Supports the sequential application of functions: \code{\link[phenex]{modelNDVI}} (used to
#' model, smooth andinterpolate VI time-series) and \code{\link[phenex]{phenoPhase}} used for
#' Phenological Phase Extraction.
#'
#'
#' @param rst A \code{RasterStack} or \code{RasterBrick} object
#'
#' @param index An integer index with the same value for each image in a year; e.g., considering three
#' images per year and a total of two years we got: c(1,1,1, 2,2,2).
#'
#' @param years An integer vector detailing which years are in the time-series. Note that the number of
#' unique values in \code{index} must be equal to the length of \code{years}.
#'
#' @param jdays Julian days for each image considering a single year (e.g., for MODIS MOD13Q1 NDVI 16-day
#' product, the following could be used \code{jdays = seq(1,365,16)}, default value).
#'
#' @param viMultFactor A multiplication factor for vegetation indices for putting data in the decimal
#' format required by package phenex.
#'
#' @param export How to export data. Either a string with \code{'memory'} which will put every single output
#' including VI modelled values and \code{phenoPhase} metrics by year in memory or, \code{'file'} which export
#' everything to GeoTIFF files. The first option (\code{'memory'}) is not suited for large datasets.
#'
#' @param methodPheno Determines whether a global or local threshold is used for greenup and senescence
#' extraction. “global” threshold: The day of the year is returned, where NDVI values are first equal or
#' higher as the value of ‘threshold’. If the threshold is higher than the values of the timeseries, ‘-1’ will
#' be returned.“local” threshold: The day of the year is returned, for which NDVI values first reach the value
#' of ‘threshold’ (interpreted as percentage) between lowest and highest NDVI value of the time-series. The lowest
#' NDVI value is chosen depending on phase selected. For “greenup”, the lowest value before day of maximum NDVI
#' value is used. For “senescence”, the lowest value after day of maximum NDVI value is used.
#'
#' @param ... Further arguments to be passed to \code{\link[phenex]{modelNDVI}}. These are parameters passed
#' to smoothing or modelling functions. These are: ‘slidingperiod’ for correction “bise”, ‘window.ravg’ for
#' correction “ravg”, ‘asym’ for method “Gauss”, ‘filter.threshold’ for method “FFT” and ‘degree’, ‘window.sav’
#' and ‘smoothing’ for method “SavGol”.
#'
#' @inherit phenex::modelNDVI params
#'
#' @inherit phenex::phenoPhase params
#'
#' @return If \code{export = 'memory'}, then a list containing the following elements will be returned:
#'
#' \describe{
#'
#'     \item{modVIRstStack}{A RasterStack containing modelled data for all years}
#'
#'     \item{phenoMetricsRstAvg}{A RasterStack containing \code{phenoPhase} metrics average}
#'
#'     \item{phenoMetricsRstStd}{A RasterStack containing \code{phenoPhase} metrics standard-deviation}
#' }
#'
#' Otherwise \code{NULL} will be returned.
#'
#' @details
#' The function supports multiple \code{phase} values (character vector) to be passed to \code{phenoPhase} function.
#'
#' @seealso
#' \code{\link[phenex:modelNDVI]{modelNDVI}}
#' \code{\link[phenex:phenoPhase]{phenoPhase}}

#'
#' @note
#' These functions are experimental and errors may happen. In addition, the performance of \code{phenoPhase} may
#' not be suitable to work with large datasets (npixels > 5E5 or so). A couple of limitations apply to this function:
#'
#' \itemize{
#'    \item Only complete years can be used (at least for now);
#'
#'    \item only a single value can be passed to \code{modelNDVI} function: \code{multipleSeasons},
#'    \code{correction}, and, \code{method} parameters; as well as for \code{phenoPhase} function: \code{method}
#'    and \code{threshold} parameters.
#' }
#'
#' @importFrom raster values<-
#' @importFrom raster values
#' @importFrom raster ncell
#' @importFrom raster writeRaster
#' @importFrom raster stack
#' @importFrom raster nlayers
#' @importFrom raster raster
#' @importFrom phenex modelNDVI
#' @importFrom phenex phenoPhase
#' @importFrom phenex leapYears
#'
#' @export

phenexRaster <- function(rst, index, years, jdays = seq(1,365,16),
                       multipleSeasons = FALSE, correction = "bise",
                       method = "LinIP", doParallel = FALSE, viMultFactor = NULL,
                       phase, methodPheno = "local", threshold = 0.55, n = 1000,
                       export = "memory", ...){

  # Check input data
  if(!inherits(rst,c("RasterStack","RasterBrick")))
    stop("rst must be an object of class RasterStack or RasterBrick")

  # Check the indices
  if(length(years)!=length(unique(index)))
    stop("Number of years must be equal to the number of unique values in index!")

  # Check options in phse
  if(!all((phase %in% c("max","maxval","min","minval","greenup","senescence"))))
    stop("phase must be one of the following options: max, maxval, min, minval, greenup or senescence")


  # Get useful metadata
  layerIndices <- 1:(raster::nlayers(rst))


  # Hold export data in lists
  if(export == "memory"){
    phenoMetricsRstAvg <- list()
    phenoMetricsRstStd <- list()
  }


  j <- 0
  for(i in unique(index)){

    j <- j + 1
    yr <- years[j]

    cat("### ----- Processing year [",yr,"] ----- ##\n\n")

    indexBool <- index == i
    nImgsYear <- sum(indexBool)

    if(nImgsYear != length(jdays))
      stop("Number of images within each year must be equal to the length of jdays
           (the Julian day of each image)!")

    nc <- ifelse(phenex::leapYears(yr), 366, 365)
    layersToProc <- layerIndices[indexBool]

    # Extract values from raster layers ----------------------------------------------------
    #
    rstValues <- raster::values(rst[[layersToProc]])

    # Multiply raster values?
    if(!is.null(viMultFactor)) rstValues <- rstValues * viMultFactor

    # Fill year matrix
    ndviMat <- matrix(NA, ncol = nc, nrow = raster::ncell(rst))

    # Make a 'pseudo-daily' NDVI matrix from composite values
    ndviMat[,jdays] <- rstValues

    # Model, smooth and interpolate VI values ----------------------------------------------
    #
    ndviObj <- phenex::modelNDVI(ndviMat,
                         year.int = yr,
                         multipleSeasons = multipleSeasons,
                         correction = correction,
                         method = method,
                         MARGIN = 1, # Apply functions over rows
                         silent = TRUE,
                         doParallel = doParallel, ...)


    # Make a matrix object from the modelled VI data
    ndviModData <- getModelledValuesToMatrix(ndviObj)

    # Subset only to julian days of the original data
    ndviModData <- ndviModData[,jdays]

    # Create a new RasterStack
    newTmpRstStack <- rst
    values(newTmpRstStack) <- ndviModData

    if(export == "memory"){

      if(j==1){
        modVIRstStack <- newTmpRstStack
      }else{
        modVIRstStack <- raster::stack(modVIRstStack, newTmpRstStack)
      }

    }else if (export == "none"){

      modVIRstStack <- NULL

    }else if(export == "file"){

      raster::writeRaster(newTmpRstStack, filename = paste("viModelledValues_yr",yr,".tif",sep=""))

    }else{
      stop("export options must be equal to: \'memory\', \'file\', or \'none\'")
    }

    # Calculate phase values ----------------------------------------------------------------
    #
    #

    for(phaseValue in phase){

      cat("Calculating phase metric:", phaseValue,"....")

      # Calculate and extract phenology metrics
      phaseValues <- lapply(ndviObj, FUN = function(x) phenex::phenoPhase(x, phaseValue, methodPheno, threshold, n))

      # Get phase values for the average and std
      phaseValuesAvg <- unlist(lapply(phaseValues, function(x) x$mean), use.names = FALSE)
      phaseValuesStd <- unlist(lapply(phaseValues, function(x) x$sd), use.names = FALSE)

      # Create new raster dataset
      phaseRstTmpAvg <- raster::raster(rst[[1]])
      phaseRstTmpStd <- raster::raster(rst[[1]])

      # Set values for selected phenometric average and std-deviation
      values(phaseRstTmpAvg) <- phaseValuesAvg
      values(phaseRstTmpStd) <- phaseValuesStd

      if(export == "memory"){

        # Accumulate rasters
        if(j==1){
          phenoMetricsRstAvg[[phaseValue]] <- phaseRstTmpAvg
          phenoMetricsRstStd[[phaseValue]] <- phaseRstTmpStd
        }else{
          phenoMetricsRstAvg[[phaseValue]] <- raster::stack(phenoMetricsRstAvg[[phaseValue]], phaseRstTmpAvg)
          phenoMetricsRstStd[[phaseValue]] <- raster::stack(phenoMetricsRstStd[[phaseValue]], phaseRstTmpStd)
        }

      }else if(export == "file"){

        raster::writeRaster(phaseRstTmpAvg, filename = paste("phenoPhaseAvg_",yr,"_",phaseValue,".tif",sep=""))
        raster::writeRaster(phaseRstTmpStd, filename = paste("phenoPhaseStd_",yr,"_",phaseValue,".tif",sep=""))


      }else{
        stop("export options must be equal to: \'memory\', \'file\', or \'none\'")
      }

      cat("done.\n")

    }

    cat("\nFinished year! \n\n")

  }


  if(export == "file"){
    return(NULL)
  }else{
    # Export final object
    return(list(

      modVIRstStack = modVIRstStack,
      phenoMetricsRstAvg = phenoMetricsRstAvg,
      phenoMetricsRstStd = phenoMetricsRstStd
    ))

  }
}