R/pvts_cd.R

In pvts-approach/PVts: photosynthetic vegetation time series manipulation in R.

#' Hamunyela smoothing
#'
#' Perform the Hamunyela smoothing (see notes) in a time series vector
#' to eliminate negative outliers.
#'
#' @note Hamunyela, E., Verbesselt, J., Roerink, G., & Herold, M. (2013).
#' Trends in spring phenology of western European deciduous forests. Remote Sensing,
#'  5(12), 6159-6179.
#'
#' @export
#' @param x Numeric vector
#' @importFrom forecast na.interp
#' @importFrom zoo rollapply
#' @examples
#' library(PVts)
#'
#' x <- c(1,2,4,-100,2,5,10)
#' pvts_hamunyela(x)
#'
pvts_hamunyela <- function(x) {
  np <- try(na.interp(x))
  if (class(np)=='try-error') {
    np <- rep(100,length(x))
  } else {
    np <- .Call(`_PVts_rcpp_hamunyela`, np)
  }
  return(np)
}


#' Estimate the mean and sd after applied the Hamunyela smoothing in a time series vector (TODO)
#'
#' @param x Numeric vector.
#'
pvts_hamunyela2 <- function(x) {
  serie_data <- pvts_hamunyela(x)
  np <- c(mean(serie_data),sd(serie_data))
  return(np)
}

#' Just estimate the mean and sd
#'
#' @param x Numeric vector.
#'
pvts_msdpixel <- function(x) {
  np <- try(na.interp(x))
  if (class(np)=='try-error') {
    np <- c(100,0)
  } else {
    stats <- c(mean(np),sd(np))
  }
  return(stats)
}


#' Mean and Standard Deviation in a Collection of Images
#'
#' Estimate the mean and the standard deviation in a Multitemporal Images.
#'
#' @export
#' @param x filename(character), Raster*  or star object (see \link[stars]{read_stars}).
#' @param CLUSTER cluster to use for parallel apply; see \link[parallel]{makeCluster}.
#' @param hamunyela logical; if \code{TRUE}, the noise will be removed applying the Hamunyela smoothing.
#' @param fit_negative logical; if \code{TRUE}, negative values are replace by NA.
#' @param RasterLayer logical; if \code{TRUE}, return a RasterLayer object.
#' @importFrom methods as is
#' @importFrom stars read_stars st_as_stars st_set_dimensions st_apply st_apply
#' @importFrom stats setNames sd
#' @importFrom raster stack
#' @examples
#' library(PVts)
#' data(pvts_datasets)
#' x = pvts_datasets$PVts[,1:50,1:50]
#'
#' f1 = pvts_meanstd(x)
#' f2 = pvts_meanstd(x,hamunyela = FALSE)
#'
pvts_meanstd <- function(x, CLUSTER = NULL, fit_negative = TRUE,
                         hamunyela = TRUE, RasterLayer = TRUE) {

  # change the object class to stars
  x <- pvts_whatkinditis(x)

  # Create a new dimension with the attributes if > 2
  if (length(dim(x)) == 2) {
    x = x %>% merge %>% st_set_dimensions(names = c("x", "y", "band"))
  }

  if (fit_negative) {
    x[x<0]=NA
  }

  if (hamunyela) {
    meansd <- x %>%
      st_apply(MARGIN = c('x','y'),
               FUN = pvts_hamunyela2,
               CLUSTER = CLUSTER) %>%
      split('pvts_hamunyela2') %>%
      setNames(c('mean','sd'))
  } else {
    meansd <- x %>%
      st_apply(MARGIN = c('x','y'),
               FUN = pvts_msdpixel,
               CLUSTER = CLUSTER) %>%
      split('pvts_msdpixel') %>%
      setNames(c('mean','sd'))
  }

  if (RasterLayer) {
    rmean <- as(meansd[1],"Raster")
    rsd <- as(meansd[2],"Raster")
    stack(rmean,rsd) %>%
      'names<-'(c('mean','sd')) -> meansd
  }

  return(meansd)
}

#' Change detection using PVts series
#'
#' Estimate the change detection considering the mean and sd and a disturbance threshold
#'
#' @export
#' @param mean filename(character), Raster*  or star object (see \link[stars]{read_stars}); represent the spatial mean.
#' @param std filename(character), Raster*  or star object (see \link[stars]{read_stars}); represent the spatial standard deviation.
#' @param Rc filename(character), Raster*  or star object (see \link[stars]{read_stars}); image to compare the change.
#' @param threshold logical; if \code{TRUE}, negative values are replace by NA.
#' @param RasterLayer logical; if \code{TRUE}, return a RasterLayer object.
#' @importFrom stars st_as_stars
#'
#' @examples
#' library(stars)
#' library(raster)
#' library(PVts)
#'
#' data(pvts_datasets)
#' x = pvts_datasets$PVts[,1:50,1:50]
#'
#' # estimate the mean and std of the first 15 years (1990-2004)
#' mean_std = pvts_meanstd(x[1:15])
#' # Get the disturbance of 2018.
#' compare = as(x[28,1:50,1:50],'Raster')
#' compare[compare<0]=NA
#' pvts_cd(mean_std[[1]],mean_std[[2]],compare)
#'
pvts_cd <- function(mean,std,Rc,threshold = 5, RasterLayer=TRUE) {

  mean <- pvts_whatkinditis(mean)
  sd <- pvts_whatkinditis(std)
  Rc <- pvts_whatkinditis(Rc)

  breakR <- Rc < (mean - threshold*sd)
  if (RasterLayer) {
    breakR <- as(breakR,"Raster")
  }
  return(breakR)
}

#' Is a character, Raster or stars object?
#'
#' @param x object to evaluate
#' @return a stars object
#'
pvts_whatkinditis <-function(x) {
  if (is.character(x)) {
      x = read_stars(x)
  } else if(is(x,'stars')) {
      x
  } else if (is(x,'RasterLayer') | is(x,'RasterStack') | is(x,'RasterBrick')) {
      x = st_as_stars(x)
  } else {
      stop(class(x), ' class is not supported')
  }
  return(x)
}