R/PLUGPhenAnoRFDPLUS.R
In MDecuy/AVOCADO: Anomaly Vegetation Change Detection
Defines functions
PLUGPhenAnoRFDMapPLUS
PLUGPhenAnoRFDPLUS
PhenRef2d
Documented in
PhenRef2d
PLUGPhenAnoRFDMapPLUS
PLUGPhenAnoRFDPLUS
#' @title Avocado algorithm - calculate likelihood of vegetation-index–time space for reference vegetation
#' @name PhenRef2d
#' @description Calculate the cumulative density distribution and their likelihood values based on reference vegetation - input for Wall-to-wall map version (i.e. PLUGPhenAnoRFDPLUS function)
#' @author Roberto O. Chavez, Mathieu Decuyper
#' @param phen Numeric vector with the values of the reference vegetation
#' @param dates A date vector. The number of dates must be equal to the number of values of time series.
#' @param h Numeric. Geographic hemisphere to define the starting date of the growing season. h=1 Northern Hemisphere; h=2 Southern Hemisphere.
#' @param anop Numeric vector with the number of values that are in x. For those values the anomalies and likelihoods will be calculated based on the phen. For example a time series has 450 values, so anop=c(1:450).
#' @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)
#' @return List which contains cumulative bivariate density distribution (cumDensity), and maximum likelihood of the vegetation-index–time space (MAXY)
#' @import npphen
#' @import stats
#' @importFrom lubridate yday
#' @seealso \code{\link{PLUGPhenAnoRFDMapPLUS}} \code{\link{PLUGPhenAnoRFDPLUS}}
#' @examples
#' \dontrun{
#' # Loading raster data
#' library(terra)
#' library(npphen)
#'
#' MDD <- rast(system.file("extdata", "MDD_NDMI_1990_2020.tif", package = "AVOCADO"))
#' # load dates vector
#' load(system.file("extdata", "MDD_dates.RData", package = "AVOCADO"))
#' # load reference forest shapefile
#' MDDref <- vect(system.file("extdata", "MDDref.gpkg", package = "AVOCADO"))
#' # Create the reference curve
#' ref.ext <- ext(MDDref)
#' ref.brick <- crop(MDD, ref.ext)
#' fin <- nrow(ref.brick) * ncol(ref.brick)
#' phen <- as.numeric(terra::extract(ref.brick, 1))
#' d1 <- MDD_dates
#' for (i in 2:fin) {
#' pp <- as.numeric(terra::extract(ref.brick, i))
#' phen <- c(phen, pp)
#' d1 <- c(d1, MDD_dates)
#' }
#' MDD_fref <- PhenRef2d(phen, d1, h = 1, anop = c(1:1063), rge = c(0, 10000))
#'
#' # plot reference curve + probabilities (same result as using PhenKplot())
#' image(MDD_fref$x, MDD_fref$y, MDD_fref$cumDensity, xlab = "DOY", ylab = "NMDI",
#' font.lab = 2, breaks = c(0, 0.5, 0.75, 0.9, 0.95), col = grDevices::heat.colors(n = 4, alpha = 0.6))
#' contour(MDD_fref$x, MDD_fref$y, MDD_fref$cumDensity, levels = c(0, 0.5, 0.75, 0.9, 0.95), add = T, col = grDevices::grey(0.25), labcex = 1)
#' lines(seq(1, 365), MDD_fref$MAXY, lwd = 3, col = "dark red")
#'
#' }
#'
#' @export
PhenRef2d <-
function(phen, dates, h, anop, rge) {
# 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(phen)) {
stop("N of dates and files do not match")
}
ano.min <- min(anop)
ano.max <- max(anop)
ano.len <- ano.max - ano.min + 1
len2 <- 2 * ano.len
if (ano.min >= ano.max) {
stop("for anop, lower value > upper value")
}
DOY <- lubridate::yday(dates)
DOY[which(DOY == 366)] <- 365
D1 <- cbind(DOY, phen)
if (length(unique(D1[, 2])) < 10 | (nrow(D1) - sum(is.na(D1))) < (0.1 * length(D1))) {
return(rep(NA, len2))
}
# 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 <- ks::Hpi(na.omit(D1))
Hmat[1, 2] <- Hmat[2, 1]
K1 <- ks::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) {
Kdiv <- sum(K1$estimate[j, ])
ifelse(Kdiv == 0, K1Con[j, ] <- 0, K1Con[j, ] <- K1$estimate[j, ] / sum(K1$estimate[j, ]))
}
first.no.NA.DOY <- min(D1[, 1][which(is.na(D1[, 2]) == FALSE)])
last.no.NA.DOY <- max(D1[, 1][which(is.na(D1[, 2]) == FALSE)])
MAXY <- apply(K1Con, 1, max)
for (i in 1:365) {
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])
}
if (i < first.no.NA.DOY) {
MAXY[i] <- NA
}
if (i > last.no.NA.DOY) {
MAXY[i] <- NA
}
}
# c. Calculating cumulative bivariate density distribution
h2d <- list()
h2d$x <- seq(1, 365)
h2d$y <- seq(rge[1], rge[2], len = 500)
h2d$density <- K1Con / sum(K1Con)
uniqueVals <- rev(unique(sort(h2d$density)))
cumRFDs <- cumsum(uniqueVals)
names(cumRFDs) <- uniqueVals
h2d$cumDensity <- matrix(nrow = nrow(h2d$density), ncol = ncol(h2d$density))
h2d$cumDensity[] <- cumRFDs[as.character(h2d$density)]
na.sta <- first.no.NA.DOY - 1
na.end <- last.no.NA.DOY + 1
if (na.sta >= 1) {
h2d$cumDensity[1:na.sta, ] <- NA
}
if (na.end <= 365) {
h2d$cumDensity[na.end:365, ] <- NA
}
# d. prepare output: attach MAXY to h2d list
h2d$MAXY <- MAXY
h2d
}
#' @title Avocado algorithm - Single pixel version
#' @name PLUGPhenAnoRFDPLUS
#' @description Calculate the Anomaly and their likelihood values (based on the difference between the pixel values and reference vegetation) for a numeric vector
#' @author Roberto O. Chavez, Mathieu Decuyper
#' @param x numeric vector. Time series of vegetation index (e.g. NDVI, EVI, NDMI)
#' @param phenref List which contains cumulative bivariate density distribution and maximum likelihood of the vegetation-index–time space based on reference vegetation obtained from \code{\link{PhenRef2d}}
#' @param dates A date vector. The number of dates must be equal to the number of values of time series.
#' @param h Numeric. Geographic hemisphere to define the starting date of the growing season. h=1 Northern Hemisphere; h=2 Southern Hemisphere.
#' @param anop Numeric vector with the number of values that are in x. For those values the anomalies and likelihoods will be calculated based on the phen. For example a time series has 450 values, so anop=c(1:450).
#' @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)
#' @return vector of length([x]*2) containing all anomalies, followed by their position within the reference frequency distribution (RFD)
#' @seealso \code{\link{PLUGPhenAnoRFDMapPLUS}}
#' @import npphen
#' @import stats
#' @importFrom lubridate yday
#' @importFrom graphics abline
#' @examples
#' \dontrun{
#' # ===================================================================================================================
#' # Loading raster data
#' library(terra)
#' library(npphen)
#' MDD <- rast(system.file("extdata", "MDD_NDMI_1990_2020.tif", package = "AVOCADO"))
#' # load dates vector
#' load(system.file("extdata", "MDD_dates.RData", package = "AVOCADO"))
#' # load reference forest data (output from PhenRef2d)
#' load(system.file("extdata", "MDD_forestReference.RData", package = "AVOCADO"))
#' ## time series extraction for a single pixel
#' px <- vect(cbind(-69.265, -12.48))
#' plot(MDD[[1]])
#' plot(px, add = T)
#'
#' # extract series
#' px_series <- as.numeric(terra::extract(MDD, px, ID = F))
#' plot(MDD_dates, px_series, type = "b", xlab = "", ylab = "NDMI")
#' ## Anomaly calculation
#' anom_rfd <- PLUGPhenAnoRFDPLUS(x = px_series, phenref = MDD_fref, dates = MDD_dates, h = 2, anop = c(1:1063), rge = c(1, 10000))
#' }
#'
#' @export
PLUGPhenAnoRFDPLUS <-
function(x, phenref, dates, h, anop, rge) {
# 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")
}
# ref.min <- min(refp)
# ref.max <- max(refp)
ano.min <- min(anop)
ano.max <- max(anop)
ano.len <- ano.max - ano.min + 1
len2 <- 2 * ano.len
if (ano.min >= ano.max) {
stop("for anop, lower value > upper value")
}
if (all(is.na(x))) {
return(rep(NA, len2))
}
DOY <- lubridate::yday(dates)
DOY[which(DOY == 366)] <- 365
D2 <- cbind(DOY[ano.min:ano.max], x[ano.min:ano.max])
if (all(is.na(D2[, 2]))) {
return(rep(NA, len2))
}
# b. Calculating anomalies AND their RFD based on D2
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(D2)) {
D2[i, 1] <- DOGS[which(DOGS[, 1] == D2[i, 1], arr.ind = TRUE), 2]
}
}
# b.1 Anomalies
MAXY <- phenref$MAXY
Anoma <- rep(NA, ano.len)
for (i in 1:nrow(D2)) {
Anoma[i] <- as.integer(D2[i, 2] - MAXY[D2[i, 1]])
}
Anoma[1:ano.len]
# b.2 RFD
rowAnom <- matrix(NA, nrow = nrow(D2), ncol = 500)
for (i in 1:nrow(D2)) {
rowAnom[i, ] <- abs(phenref$y - D2[i, 2])
}
rowAnom2 <- unlist(apply(rowAnom, 1, function(x) {
if (all(is.na(x))) {
NA
} else {
which.min(x)
}
}))
AnomRFD <- rep(NA, nrow(D2))
for (i in 1:nrow(D2)) {
AnomRFD[i] <- phenref$cumDensity[D2[i, 1], rowAnom2[i]]
}
AnomRFDPerc <- round(100 * (AnomRFD))
plot(AnomRFDPerc[1:ano.len], xlab = "Time", ylab = "RFD (%)", font.lab = 2)
abline(h = 90, col = "red")
abline(h = 95, col = "red")
abline(h = 99, col = "red")
AnomRFDPerc[1:ano.len]
# Two results in a single numeric vector
AnomPlusRFD <- c(Anoma[1:ano.len], AnomRFDPerc[1:ano.len])
AnomPlusRFD[1:len2]
}
#' @title Avocado algorithm - Wall-to-wall map version
#' @name PLUGPhenAnoRFDMapPLUS
#' @description Calculate the Anomaly and their likelihood values (based on the difference between the pixel values and reference vegetation) for a RasterStack
#' @author Roberto O. Chavez, Mathieu Decuyper
#' @param s SpatRaster time series of vegetation index (e.g. NDVI, EVI, NDMI)
#' @param anop Numeric vector with the number of values that are in x. For those values the anomalies and likelihoods will be calculated based on the phen. For example a time series has 450 values, so anop=c(1:450).
#' @param phenref List which contains cumulative bivariate density distribution and maximum likelihood of the vegetation-index–time space based on reference vegetation obtained from \code{\link{PhenRef2d}}
#' @param dates A date vector. The number of dates must be equal to the number of values of time series.
#' @param h Numeric. Geographic hemisphere to define the starting date of the growing season. h=1 Northern Hemisphere; h=2 Southern Hemisphere.
#' @param nCluster Numeric. Number of CPU's to be used for the job.
#' @param outname Character vector with the output path and filename with extension or only the filename and extension if work directory was set. More information: See writeRaster
#' @param datatype Character vector that determines the interpretation of values written to disk. More information: See \code{\link[terra]{writeRaster}}
#' @return SpatRaster of nlayers([s]stack*2) containing all anomalies, followed by their likelihoods
#' @import npphen
#' @import terra
#' @import stats
#' @importFrom lubridate yday
#' @seealso \code{\link{PLUGPhenAnoRFDPLUS}}
#' @examples
#' \dontrun{
#' # Loading raster data
#' library(terra)
#' library(npphen)
#'
#' MDD <- rast(system.file("extdata", "MDD_NDMI_1990_2020.tif", package = "AVOCADO"))
#' # load dates vector
#' load(system.file("extdata", "MDD_dates.RData", package = "AVOCADO"))
#' # load reference forest data (output from PhenRef2d)
#' load(system.file("extdata", "MDD_forestReference.RData", package = "AVOCADO"))
#'
#' ## Anomaly calculation
#' # checking availiable cores and leave one free
#' nc1 <- parallel::detectCores() - 1
#' PLUGPhenAnoRFDMapPLUS(
#' s = MDD, dates = MDD_dates, h = 1, phenref = MDD_fref, anop = c(1:1063),
#' nCluster = nc1, outname = "YourDirectory/MDD_AnomalyLikelihood.tif",
#' datatype = "INT2S")
#' )
#' # The output file contains all the anomalies, followed by all their likelihoods
#' # and thus has twice the number of layers as the input raster stack.
#' }
#' @export
PLUGPhenAnoRFDMapPLUS <-
function(s, phenref, dates, h, anop, nCluster, outname, datatype) {
ff <- function(x) {
# a.Preparing dataset
if (length(dates) != length(x)) {
stop("N of dates and files do not match")
}
if (length(x) < length(anop)) {
stop("Inconsistent anop. Argument anop can't be grater than length(x)")
}
# ref.min <- min(refp)
# ref.max <- max(refp)
ano.min <- min(anop)
ano.max <- max(anop)
ano.len <- ano.max - ano.min + 1
len2 <- 2 * ano.len
if (ano.min >= ano.max) {
stop("for anop, lower value > upper value")
}
if (all(is.na(x))) {
return(rep(NA, len2))
}
DOY <- lubridate::yday(dates)
DOY[which(DOY == 366)] <- 365
D2 <- cbind(DOY[ano.min:ano.max], x[ano.min:ano.max])
if (all(is.na(D2[, 2]))) {
return(rep(NA, len2))
}
if (h != 1 && h != 2) {
stop("Invalid h")
}
DOGS <- cbind(seq(1, 365), c(seq(185, 365), seq(1, 184)))
# b. Calculating anomalies AND their RFD based on D2
if (h == 2) {
for (i in 1:nrow(D2)) {
D2[i, 1] <- DOGS[which(DOGS[, 1] == D2[i, 1], arr.ind = TRUE), 2]
}
}
# b.1 Anomalies
MAXY <- phenref$MAXY
Anoma <- rep(NA, ano.len)
for (i in 1:nrow(D2)) {
Anoma[i] <- as.integer(D2[i, 2] - MAXY[D2[i, 1]])
}
Anoma[1:ano.len]
# b.2 RFD
rowAnom <- matrix(NA, nrow = nrow(D2), ncol = 500)
for (i in 1:nrow(D2)) {
rowAnom[i, ] <- abs(phenref$y - D2[i, 2])
}
rowAnom2 <- unlist(apply(rowAnom, 1, function(x) {
if (all(is.na(x))) {
NA
} else {
which.min(x)
}
}))
AnomRFD <- rep(NA, nrow(D2))
for (i in 1:nrow(D2)) {
AnomRFD[i] <- phenref$cumDensity[D2[i, 1], rowAnom2[i]]
}
AnomRFDPerc <- round(100 * (AnomRFD))
AnomRFDPerc[1:ano.len]
# Two results in a single numeric vector
AnomPLUSRFD <- c(Anoma[1:ano.len], AnomRFDPerc[1:ano.len])
AnomPLUSRFD[1:len2]
}
#----------------------------------------------------------------------------------------
# cluster processing
app(s, fun = ff, filename = outname, cores = nCluster, overwrite = T, wopt = list(datatype = datatype))
}
MDecuy/AVOCADO documentation built on April 14, 2025, 5:30 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions PLUGPhenAnoRFDMapPLUS PLUGPhenAnoRFDPLUS PhenRef2d
Documented in PhenRef2d PLUGPhenAnoRFDMapPLUS PLUGPhenAnoRFDPLUS
#' @title Avocado algorithm - calculate likelihood of vegetation-index–time space for reference vegetation
#' @name PhenRef2d
#' @description Calculate the cumulative density distribution and their likelihood values based on reference vegetation - input for Wall-to-wall map version (i.e. PLUGPhenAnoRFDPLUS function)
#' @author Roberto O. Chavez, Mathieu Decuyper
#' @param phen Numeric vector with the values of the reference vegetation
#' @param dates A date vector. The number of dates must be equal to the number of values of time series.
#' @param h Numeric. Geographic hemisphere to define the starting date of the growing season. h=1 Northern Hemisphere; h=2 Southern Hemisphere.
#' @param anop Numeric vector with the number of values that are in x. For those values the anomalies and likelihoods will be calculated based on the phen. For example a time series has 450 values, so anop=c(1:450).
#' @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)
#' @return List which contains cumulative bivariate density distribution (cumDensity), and maximum likelihood of the vegetation-index–time space (MAXY)
#' @import npphen
#' @import stats
#' @importFrom lubridate yday
#' @seealso \code{\link{PLUGPhenAnoRFDMapPLUS}} \code{\link{PLUGPhenAnoRFDPLUS}}
#' @examples
#' \dontrun{
#' # Loading raster data
#' library(terra)
#' library(npphen)
#'
#' MDD <- rast(system.file("extdata", "MDD_NDMI_1990_2020.tif", package = "AVOCADO"))
#' # load dates vector
#' load(system.file("extdata", "MDD_dates.RData", package = "AVOCADO"))
#' # load reference forest shapefile
#' MDDref <- vect(system.file("extdata", "MDDref.gpkg", package = "AVOCADO"))
#' # Create the reference curve
#' ref.ext <- ext(MDDref)
#' ref.brick <- crop(MDD, ref.ext)
#' fin <- nrow(ref.brick) * ncol(ref.brick)
#' phen <- as.numeric(terra::extract(ref.brick, 1))
#' d1 <- MDD_dates
#' for (i in 2:fin) {
#' pp <- as.numeric(terra::extract(ref.brick, i))
#' phen <- c(phen, pp)
#' d1 <- c(d1, MDD_dates)
#' }
#' MDD_fref <- PhenRef2d(phen, d1, h = 1, anop = c(1:1063), rge = c(0, 10000))
#'
#' # plot reference curve + probabilities (same result as using PhenKplot())
#' image(MDD_fref$x, MDD_fref$y, MDD_fref$cumDensity, xlab = "DOY", ylab = "NMDI",
#' font.lab = 2, breaks = c(0, 0.5, 0.75, 0.9, 0.95), col = grDevices::heat.colors(n = 4, alpha = 0.6))
#' contour(MDD_fref$x, MDD_fref$y, MDD_fref$cumDensity, levels = c(0, 0.5, 0.75, 0.9, 0.95), add = T, col = grDevices::grey(0.25), labcex = 1)
#' lines(seq(1, 365), MDD_fref$MAXY, lwd = 3, col = "dark red")
#'
#' }
#'
#' @export
PhenRef2d <-
function(phen, dates, h, anop, rge) {
# 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(phen)) {
stop("N of dates and files do not match")
}
ano.min <- min(anop)
ano.max <- max(anop)
ano.len <- ano.max - ano.min + 1
len2 <- 2 * ano.len
if (ano.min >= ano.max) {
stop("for anop, lower value > upper value")
}
DOY <- lubridate::yday(dates)
DOY[which(DOY == 366)] <- 365
D1 <- cbind(DOY, phen)
if (length(unique(D1[, 2])) < 10 | (nrow(D1) - sum(is.na(D1))) < (0.1 * length(D1))) {
return(rep(NA, len2))
}
# 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 <- ks::Hpi(na.omit(D1))
Hmat[1, 2] <- Hmat[2, 1]
K1 <- ks::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) {
Kdiv <- sum(K1$estimate[j, ])
ifelse(Kdiv == 0, K1Con[j, ] <- 0, K1Con[j, ] <- K1$estimate[j, ] / sum(K1$estimate[j, ]))
}
first.no.NA.DOY <- min(D1[, 1][which(is.na(D1[, 2]) == FALSE)])
last.no.NA.DOY <- max(D1[, 1][which(is.na(D1[, 2]) == FALSE)])
MAXY <- apply(K1Con, 1, max)
for (i in 1:365) {
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])
}
if (i < first.no.NA.DOY) {
MAXY[i] <- NA
}
if (i > last.no.NA.DOY) {
MAXY[i] <- NA
}
}
# c. Calculating cumulative bivariate density distribution
h2d <- list()
h2d$x <- seq(1, 365)
h2d$y <- seq(rge[1], rge[2], len = 500)
h2d$density <- K1Con / sum(K1Con)
uniqueVals <- rev(unique(sort(h2d$density)))
cumRFDs <- cumsum(uniqueVals)
names(cumRFDs) <- uniqueVals
h2d$cumDensity <- matrix(nrow = nrow(h2d$density), ncol = ncol(h2d$density))
h2d$cumDensity[] <- cumRFDs[as.character(h2d$density)]
na.sta <- first.no.NA.DOY - 1
na.end <- last.no.NA.DOY + 1
if (na.sta >= 1) {
h2d$cumDensity[1:na.sta, ] <- NA
}
if (na.end <= 365) {
h2d$cumDensity[na.end:365, ] <- NA
}
# d. prepare output: attach MAXY to h2d list
h2d$MAXY <- MAXY
h2d
}
#' @title Avocado algorithm - Single pixel version
#' @name PLUGPhenAnoRFDPLUS
#' @description Calculate the Anomaly and their likelihood values (based on the difference between the pixel values and reference vegetation) for a numeric vector
#' @author Roberto O. Chavez, Mathieu Decuyper
#' @param x numeric vector. Time series of vegetation index (e.g. NDVI, EVI, NDMI)
#' @param phenref List which contains cumulative bivariate density distribution and maximum likelihood of the vegetation-index–time space based on reference vegetation obtained from \code{\link{PhenRef2d}}
#' @param dates A date vector. The number of dates must be equal to the number of values of time series.
#' @param h Numeric. Geographic hemisphere to define the starting date of the growing season. h=1 Northern Hemisphere; h=2 Southern Hemisphere.
#' @param anop Numeric vector with the number of values that are in x. For those values the anomalies and likelihoods will be calculated based on the phen. For example a time series has 450 values, so anop=c(1:450).
#' @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)
#' @return vector of length([x]*2) containing all anomalies, followed by their position within the reference frequency distribution (RFD)
#' @seealso \code{\link{PLUGPhenAnoRFDMapPLUS}}
#' @import npphen
#' @import stats
#' @importFrom lubridate yday
#' @importFrom graphics abline
#' @examples
#' \dontrun{
#' # ===================================================================================================================
#' # Loading raster data
#' library(terra)
#' library(npphen)
#' MDD <- rast(system.file("extdata", "MDD_NDMI_1990_2020.tif", package = "AVOCADO"))
#' # load dates vector
#' load(system.file("extdata", "MDD_dates.RData", package = "AVOCADO"))
#' # load reference forest data (output from PhenRef2d)
#' load(system.file("extdata", "MDD_forestReference.RData", package = "AVOCADO"))
#' ## time series extraction for a single pixel
#' px <- vect(cbind(-69.265, -12.48))
#' plot(MDD[[1]])
#' plot(px, add = T)
#'
#' # extract series
#' px_series <- as.numeric(terra::extract(MDD, px, ID = F))
#' plot(MDD_dates, px_series, type = "b", xlab = "", ylab = "NDMI")
#' ## Anomaly calculation
#' anom_rfd <- PLUGPhenAnoRFDPLUS(x = px_series, phenref = MDD_fref, dates = MDD_dates, h = 2, anop = c(1:1063), rge = c(1, 10000))
#' }
#'
#' @export
PLUGPhenAnoRFDPLUS <-
function(x, phenref, dates, h, anop, rge) {
# 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")
}
# ref.min <- min(refp)
# ref.max <- max(refp)
ano.min <- min(anop)
ano.max <- max(anop)
ano.len <- ano.max - ano.min + 1
len2 <- 2 * ano.len
if (ano.min >= ano.max) {
stop("for anop, lower value > upper value")
}
if (all(is.na(x))) {
return(rep(NA, len2))
}
DOY <- lubridate::yday(dates)
DOY[which(DOY == 366)] <- 365
D2 <- cbind(DOY[ano.min:ano.max], x[ano.min:ano.max])
if (all(is.na(D2[, 2]))) {
return(rep(NA, len2))
}
# b. Calculating anomalies AND their RFD based on D2
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(D2)) {
D2[i, 1] <- DOGS[which(DOGS[, 1] == D2[i, 1], arr.ind = TRUE), 2]
}
}
# b.1 Anomalies
MAXY <- phenref$MAXY
Anoma <- rep(NA, ano.len)
for (i in 1:nrow(D2)) {
Anoma[i] <- as.integer(D2[i, 2] - MAXY[D2[i, 1]])
}
Anoma[1:ano.len]
# b.2 RFD
rowAnom <- matrix(NA, nrow = nrow(D2), ncol = 500)
for (i in 1:nrow(D2)) {
rowAnom[i, ] <- abs(phenref$y - D2[i, 2])
}
rowAnom2 <- unlist(apply(rowAnom, 1, function(x) {
if (all(is.na(x))) {
NA
} else {
which.min(x)
}
}))
AnomRFD <- rep(NA, nrow(D2))
for (i in 1:nrow(D2)) {
AnomRFD[i] <- phenref$cumDensity[D2[i, 1], rowAnom2[i]]
}
AnomRFDPerc <- round(100 * (AnomRFD))
plot(AnomRFDPerc[1:ano.len], xlab = "Time", ylab = "RFD (%)", font.lab = 2)
abline(h = 90, col = "red")
abline(h = 95, col = "red")
abline(h = 99, col = "red")
AnomRFDPerc[1:ano.len]
# Two results in a single numeric vector
AnomPlusRFD <- c(Anoma[1:ano.len], AnomRFDPerc[1:ano.len])
AnomPlusRFD[1:len2]
}
#' @title Avocado algorithm - Wall-to-wall map version
#' @name PLUGPhenAnoRFDMapPLUS
#' @description Calculate the Anomaly and their likelihood values (based on the difference between the pixel values and reference vegetation) for a RasterStack
#' @author Roberto O. Chavez, Mathieu Decuyper
#' @param s SpatRaster time series of vegetation index (e.g. NDVI, EVI, NDMI)
#' @param anop Numeric vector with the number of values that are in x. For those values the anomalies and likelihoods will be calculated based on the phen. For example a time series has 450 values, so anop=c(1:450).
#' @param phenref List which contains cumulative bivariate density distribution and maximum likelihood of the vegetation-index–time space based on reference vegetation obtained from \code{\link{PhenRef2d}}
#' @param dates A date vector. The number of dates must be equal to the number of values of time series.
#' @param h Numeric. Geographic hemisphere to define the starting date of the growing season. h=1 Northern Hemisphere; h=2 Southern Hemisphere.
#' @param nCluster Numeric. Number of CPU's to be used for the job.
#' @param outname Character vector with the output path and filename with extension or only the filename and extension if work directory was set. More information: See writeRaster
#' @param datatype Character vector that determines the interpretation of values written to disk. More information: See \code{\link[terra]{writeRaster}}
#' @return SpatRaster of nlayers([s]stack*2) containing all anomalies, followed by their likelihoods
#' @import npphen
#' @import terra
#' @import stats
#' @importFrom lubridate yday
#' @seealso \code{\link{PLUGPhenAnoRFDPLUS}}
#' @examples
#' \dontrun{
#' # Loading raster data
#' library(terra)
#' library(npphen)
#'
#' MDD <- rast(system.file("extdata", "MDD_NDMI_1990_2020.tif", package = "AVOCADO"))
#' # load dates vector
#' load(system.file("extdata", "MDD_dates.RData", package = "AVOCADO"))
#' # load reference forest data (output from PhenRef2d)
#' load(system.file("extdata", "MDD_forestReference.RData", package = "AVOCADO"))
#'
#' ## Anomaly calculation
#' # checking availiable cores and leave one free
#' nc1 <- parallel::detectCores() - 1
#' PLUGPhenAnoRFDMapPLUS(
#' s = MDD, dates = MDD_dates, h = 1, phenref = MDD_fref, anop = c(1:1063),
#' nCluster = nc1, outname = "YourDirectory/MDD_AnomalyLikelihood.tif",
#' datatype = "INT2S")
#' )
#' # The output file contains all the anomalies, followed by all their likelihoods
#' # and thus has twice the number of layers as the input raster stack.
#' }
#' @export
PLUGPhenAnoRFDMapPLUS <-
function(s, phenref, dates, h, anop, nCluster, outname, datatype) {
ff <- function(x) {
# a.Preparing dataset
if (length(dates) != length(x)) {
stop("N of dates and files do not match")
}
if (length(x) < length(anop)) {
stop("Inconsistent anop. Argument anop can't be grater than length(x)")
}
# ref.min <- min(refp)
# ref.max <- max(refp)
ano.min <- min(anop)
ano.max <- max(anop)
ano.len <- ano.max - ano.min + 1
len2 <- 2 * ano.len
if (ano.min >= ano.max) {
stop("for anop, lower value > upper value")
}
if (all(is.na(x))) {
return(rep(NA, len2))
}
DOY <- lubridate::yday(dates)
DOY[which(DOY == 366)] <- 365
D2 <- cbind(DOY[ano.min:ano.max], x[ano.min:ano.max])
if (all(is.na(D2[, 2]))) {
return(rep(NA, len2))
}
if (h != 1 && h != 2) {
stop("Invalid h")
}
DOGS <- cbind(seq(1, 365), c(seq(185, 365), seq(1, 184)))
# b. Calculating anomalies AND their RFD based on D2
if (h == 2) {
for (i in 1:nrow(D2)) {
D2[i, 1] <- DOGS[which(DOGS[, 1] == D2[i, 1], arr.ind = TRUE), 2]
}
}
# b.1 Anomalies
MAXY <- phenref$MAXY
Anoma <- rep(NA, ano.len)
for (i in 1:nrow(D2)) {
Anoma[i] <- as.integer(D2[i, 2] - MAXY[D2[i, 1]])
}
Anoma[1:ano.len]
# b.2 RFD
rowAnom <- matrix(NA, nrow = nrow(D2), ncol = 500)
for (i in 1:nrow(D2)) {
rowAnom[i, ] <- abs(phenref$y - D2[i, 2])
}
rowAnom2 <- unlist(apply(rowAnom, 1, function(x) {
if (all(is.na(x))) {
NA
} else {
which.min(x)
}
}))
AnomRFD <- rep(NA, nrow(D2))
for (i in 1:nrow(D2)) {
AnomRFD[i] <- phenref$cumDensity[D2[i, 1], rowAnom2[i]]
}
AnomRFDPerc <- round(100 * (AnomRFD))
AnomRFDPerc[1:ano.len]
# Two results in a single numeric vector
AnomPLUSRFD <- c(Anoma[1:ano.len], AnomRFDPerc[1:ano.len])
AnomPLUSRFD[1:len2]
}
#----------------------------------------------------------------------------------------
# cluster processing
app(s, fun = ff, filename = outname, cores = nCluster, overwrite = T, wopt = list(datatype = datatype))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.