R/PhenMap.R

In JoseLastra/npphen: Vegetation Phenological Cycle and Anomaly Detection using Remote Sensing Data

#' @title PhenMap
#' @description Estimates annual Land Surface Phenology (LSP) using time series of  a vegetation greenness raster stack.
#' @encoding UTF-8
#' @param s Raster stack with greenness (e,g. NDVI or EV) values
#' @param dates Vector with dates at which the greenness values were recorded
#' @param h Numeric. Indicates the geographic hemisphere to define the starting date of the growing season. h=1 if the vegetation is in the Northern Hemisphere (season starting at January 1st), h=2 if it is in the Southern Hemisphere (season starting at July 1st)
#' @param frequency Character string. Defines the number of samples for the output phenology and must be one of the this: *'daily'* (default) giving output vector of length 365, *'8-days'* giving output vector of length 46 (i.e MOD13Q1 and MYD13Q1), *'monthly'* giving output vector of length 12,*'bi-weekly'* giving output vector of length 24 (i.e. GIMMS) or *'16-days'* giving output vector of length 23 (i.e MOD13Q1 or MYD13Q1).
#' @param nCluster Numeric. Number of CPU cores to be used for computational calculations
#' @param outname Character vector with the output path and filename with extension or only the filename and extension if work directory was set. For example outname="output_phen.tif". See \code{\link{writeRaster}}
#' @param format Character. Output file type. See \code{\link{writeFormats}}
#' @param datatype Character. Output data type. See \code{\link{dataType}}
#' @param rge  A vector containing minimum and maximum values of the response variable used in the analysis. We suggest the use of theoretically based limits. For example in the case of MODIS NDVI or EVI, it ranges from 0 to 10,000, so rge =c(0,10000)
#' @details Derives the annual Land Surface Phenological (LSP) cycle for a standard growing season using a raster stack of satellite based greenness values such as the Normalized Difference Vegetation Index (NDVI) or Enhanced Vegetation Index (EVI). The LSP cycle is calculated for all pixels of the input raster stack in the same way as for the Phen function. The output is a multiband raster where every band is the expected greeness value at a given time step of the standard growing season. For example, for MODIS Vegetation Index 16-days composites the number of time steps of the growing season (nGS) is 23 , and therefore, the output raster will have 23 bands. A vector with dates for the greenness values is also required.
#' @return RasterStack
#' @seealso \code{\link{Phen}}
#' @examples
#' \dontrun{
#' ##DEPENDING ON HARDWARE, THIS PROCESS CAN BE HIGHLY TIME CONSUMING##
#'
#' ## Testing raster data from Central Chile (NDVI), h=2##
#'
#' # Load data
#' #RasterStack
#' data("MegaDrought_stack")
#' #Dates
#' data("dates")
#'
#' # Making the LSP raster, n bands = 23
#'
#' library(snow)
#'
#' # Define the number of cores to be use. In this example we use 1
#' nc1<-1
#'
#'PhenMap(s = MegaDrought_stack,dates = dates,h = 2,
#'frequency = '16-days',nCluster = nc1,outname = 'phen_MD.tif',
#'format = 'GTiff',datatype = 'INT2S',rge = c(0,1000))
#' #map1 <- stack("phen_MD.tif")#run only for load phenology brick
#' #plot(map1)
#'
#' ## Testing raster data from Blooming desert, Northern Chile (NDVI), h=2 ##
#'
#' # Load data
#' RasterStack
#' data("Bdesert_stack")
#' #Dates
#' data("dates")
#'
#' # Making the LSP raster, n bands = 23
#' # Define the number of cores to be use. In this example we use 1
#' nc1<-1
#'
#' PhenMap(s= Bdesert_stack,dates=dates,h=2,
#' frequency='16-days', nCluster=nc1,outname="phen_BD.tif", 
#' format="GTiff", datatype="INT2S",rge=c(0,10000))
#' #map2 <- stack("phen_BD.tif") #run only for load phenology brick
#' #plot(map2)
#' }
#' @export

PhenMap <-
  function(s,dates,h,frequency='16-days',nCluster,outname,format,datatype,rge) {
    ff <- function(x) {
      
      # a.Preparing dataset
      if(length(rge)!=2){stop("rge must be a vector of length 2")}
      if(rge[1]>rge[2]){stop("rge vector order must be minimum/maximum")}
      if(length(dates)!=length(x)){stop("N of dates and files do not match")}
      if(frequency == 'daily'){
        nGS <- 365
      }
      if(frequency == '8-days'){
        nGS <- 46
      }
      if(frequency == '16-days'){
        nGS <- 23
      }
      if(frequency == 'monthly'){
        nGS <- 12
      }
      if(frequency == 'bi-monthly'){
        nGS <- 24
      }
      if (all(is.na(x))) {
        return(rep(NA,nGS))
      }
      
      freq.method <- match(frequency, c("daily", "8-days", "16-days", 
                                         "bi-weekly", "monthly"))
      if (is.na(freq.method)){ 
        stop("Invalid frequency. Must be one of: daily, 8-days, 16-days, bi-weekly or monthly")}
      if (all(is.na(x))) {
        return(rep(NA,nGS))
      }
      
      DOY <- yday(dates)
      DOY[which(DOY==366)]<-365	#fixed leap-years
      D1<-cbind(DOY,x)
      if(length(unique(D1[,2]))<10 | (nrow(D1)-sum(is.na(D1)))<(0.1*nrow(D1))) { #replace length by nrow 
        return(rep(NA,nGS))
      }
      
      # b. Kernel calculation using the reference period (D1)
      if(h!=1 && h!=2){stop("Invalid h")}
      DOGS<-cbind(seq(1,365),c(seq(185,365),seq(1,184)))
      if(h==2){
        for(i in 1:nrow(D1)){
          D1[i,1]<-DOGS[which(DOGS[,1]==D1[i,1],arr.ind=TRUE),2]}}
      
      Hmat<-Hpi(na.omit(D1))
      Hmat[1,2]<-Hmat[2,1]
      K1<-kde(na.omit(D1),H=Hmat,xmin=c(1,rge[1]),xmax=c(365,rge[2]),gridsize=c(365,500))
      K1Con<-K1$estimate
      
      for(j in 1:365){ #Fixed problem for non calculated probs
        Kdiv<-sum(K1$estimate[j,])
        ifelse(Kdiv==0,K1Con[j,]<-0,K1Con[j,]<-K1$estimate[j,]/sum(K1$estimate[j,]))}
      
      MAXY<-apply(K1Con,1,max)
      for(i in 1:365){ #fixed artifact maxy values and blank accordingly
        n.select <- which(K1Con[i,]==MAXY[i],arr.ind=TRUE)
        if(length(n.select) > 1){
          n <- n.select[1]
          MAXY[i]<-NA}
        if(length(n.select) == 1){
          n <- n.select
          MAXY[i]<-median(K1$eval.points[[2]][n])
        }
      }
      
      # c. Getting the most probable values for a standard year (nGS time steps)
      
      if(frequency == 'daily'){
        select_DGS <- seq(1,365,1)
        Ref <- MAXY
      }
      if(frequency == '8-days'){
        select_DGS <- seq(1,365,8)
        Ref <- MAXY[select_DGS]
      }
      if(frequency == '16-days'){
        select_DGS <- seq(1,365,16)
        Ref <- MAXY[select_DGS]
      }
      if(frequency == 'monthly'){
        select_DGS <- c(15,46,74,105,135,166,196,227,258,288,319,349)
        Ref <- MAXY[select_DGS]
      }
      if(frequency == 'bi-monthly'){
        select_DGS <- c(1,15,32,46,60,74,91,105,121,135,152,166,182,196,213,227,244,258,274,288,305,319,335,349)
        Ref <- MAXY[select_DGS]
      }
      if(h==1){id.label <- 'DOY'}
      if(h==2){id.label <- 'DGS'}

      names(Ref) <- select_DGS # set names
      Ref
    }

    #----------------------------------------------------------------------------------------
    # Calculating the annual phenological curve using n clusters
    beginCluster(n=nCluster) # write 'beginCluster(n=3)' for using e.g. 3 cores, default uses all available cores)
    dates<<-dates
    clusterR(x=s,calc, args=list(ff),export=c('dates'),filename=outname,format=format,datatype=datatype,overwrite=T)
    endCluster()
}