Nothing
#' Convert a series of raster files to a single phenology raster.
#'
#' @param File_List List of raster files
#' @param PhenoFactor Character string; type of dataset to analyze (e.g., "VI", "Snow")
#' @param phase Character string; name of phenophase to be measured (e.g., "greenup", "snowmelt", "senescence" or other arguments passed to phenex::phenophase())
#' @param threshold Float threshold GWI value to be projected. Use only for VI option.
#' @param year Integer Year (YYYY)
#' @param NDVI Integer Band number of NDVI band in raster files
#' @param VIQ Integer Band number of VI Quality layer in raster files
#' @param DOY Integer Band number of Composite Day of Year layer in raster files
#' @param PR Integer Band Number of PR layer in raster files
#' @param SnowExtent Integer Band number of Maximum_Snow_Extent in raster files
#' @param verbose TRUE or FALSE (Default = FALSE)
#' @return Raster object with extent=extent(terra::rast(File_List)[1]) and CRS = crs(terra::rast(File_List)[1]). Digital numbers are expressed as Day of Year.
#' @examples
#' \dontshow{
#' fpath <- system.file("extdata", package="phenomap")
#' File_List <- paste(fpath, list.files(path = fpath, pattern=c("ExCrop_")), sep="/")
#' File_List <- File_List[1:2]
#'
#' PhenoFactor = "VI"
#' phase = "greenup"
#' threshold = 0.5
#' year = 2016
#' NDVI = 1
#' VIQ = 3
#' DOY = 4
#' PR = 5
#' verbose = TRUE
#'
#' Sample.Greenup <- mapPheno(File_List = File_List, PhenoFactor = PhenoFactor,
#' phase = phase, threshold = threshold, year = year,
#' NDVI = NDVI, VIQ = VIQ, DOY = DOY, PR = PR,
#' SnowExtent=SnowExtent,
#' verbose = verbose)
#' }
#' \dontrun{
#' fpath <- system.file("extdata", package="phenomap")
#' File_List <- paste(fpath, list.files(path = fpath, pattern=c("TinyCrop_")), sep="/")
#' File_List
#'
#' PhenoFactor = "VI"
#' phase = "greenup"
#' threshold = 0.5
#' year = 2016
#' NDVI = 1
#' VIQ = 3
#' DOY = 4
#' PR = 5
#' verbose = TRUE
#'
#' Sample.Greenup <- mapPheno(File_List = File_List, PhenoFactor = PhenoFactor,
#' phase = phase, threshold = threshold, year = year,
#' NDVI = NDVI, VIQ = VIQ, DOY = DOY, PR = PR,
#' SnowExtent=SnowExtent,
#' verbose = verbose)
#' }
#' @import stringr
#' @import terra
#' @importFrom plyr aaply
#' @importFrom phenex modelNDVI
#' @importFrom phenex phenoPhase
#' @export
mapPheno<- function(File_List = NA, PhenoFactor = NA,
phase = NA, threshold = NA, year = NA,
NDVI = NA, VIQ = NA, DOY = NA, PR = NA,
SnowExtent = NA,
verbose = FALSE){
annualcrops <- File_List
pos <- new.env()
## This function transforms a list of strings returned by searching for objects
## into a list of objects (required for proper stacking)
if(verbose){print(paste0("Sorting data into 3-dimensional arrays... ", Sys.time()))}
if(PhenoFactor == "VI"){
if(verbose){print(paste0("VI option selected... ", Sys.time()))}
if(is.na(NDVI) || is.na(VIQ) || is.na(DOY) || is.na(PR)){
stop("Crucial Dataset missing (have you included NDVI, VIQ, DOY, and PR layers?)")
}
if(!is.na(NDVI)){
if(verbose){
print(paste0("Creating NDVI array... ", Sys.time()))
print("Array dimensions... ")
print(c(dim(terra::as.array(terra::rast(annualcrops[1])))[1],
dim(terra::as.array(terra::rast(annualcrops[1])))[2],
length(annualcrops)))
}
NDVI.Array <- array(data=NA,
dim=c(dim(terra::as.array(terra::rast(annualcrops[1])))[1],
dim(terra::as.array(terra::rast(annualcrops[1])))[2],
length(annualcrops)),
dimnames=NULL)
if(verbose){print("Created blank array")}
for(i in 1:length(annualcrops)){
NDVI.Array[,,i] <- terra::as.array(terra::rast(annualcrops[i], lyrs=NDVI))
}; if(verbose){print("Filled array")}
#NDVI.Stack <- stack(annualcrops, bands = NDVI)
# if(verbose){
# print(class(NDVI.Stack))
# print(c("NDVI.Stack dimensions ", dim(NDVI.Stack)))
# }
#NDVI.Array <- as.array(NDVI.Stack)
# if(verbose){
# print("Sample NDVI matrix dimensions")
# print(dim(terra::as.array(terra::rast(annualcrops[6], band=NDVI))))
# print("NDVI.Array dimensions")
# print(dim(NDVI.Array))
# print("3")
# }
## If NDVI is less than 0, it needs to be brought up to 0
NDVIcleaner <- function(x) {x[x<0] <- 0; x};if(verbose){print("cleaner created")}
NDVIcleaned <- plyr::aaply(.data = NDVI.Array,
.margins = c(1,2,3),
.fun = NDVIcleaner)
}
if(verbose){print("NDVI cleaned")}
if(!is.na(VIQ)){
if(verbose){print(paste0("Creating VI Quality array... ", Sys.time()))}
VIQ.Array <- array(dim=c(dim(terra::as.array(terra::rast(annualcrops[1])))[1],
dim(terra::as.array(terra::rast(annualcrops[1])))[2],
length(annualcrops)),
data=NA); if(verbose){print("Created blank array")}
for(i in 1:length(annualcrops)){
VIQ.Array[,,i] <- terra::as.array(terra::rast(annualcrops[i], lyrs=VIQ))
}
if(verbose){
print("Filled VIQ Array")
print("VIQ.Array dimensions")
print(dim(VIQ.Array))
}
## first_k_bits derived from:
## https://gis.stackexchange.com/questions/144441/how-can-i-parse-modis-mod13q1-quality-layers-in-r/144487
first_k_bits <- function(int, k=16) {
## https://lpdaac.usgs.gov/products/modis_products_table/mod13q1 TABLE 2:
## MOD13Q1 VI Quality: "Bit 0 is the least significant (read bit words right to left)"
integer_vector <- as.integer(intToBits(int))[1:k]
return(paste(as.character(integer_vector),
collapse=""))
}
VIQbinary <- plyr::aaply(.data = VIQ.Array,
.margins = c(1,2,3),
.fun = first_k_bits); if(verbose){print("4")}
}
if(!is.na(DOY)){
if(verbose){print(paste0("Creating Day of Year array... ", Sys.time()))}
DOY.Array <- array(dim=c(dim(terra::as.array(terra::rast(annualcrops[1])))[1],
dim(terra::as.array(terra::rast(annualcrops[1])))[2],
length(annualcrops)),
data=NA)
for(i in 1:length(annualcrops)){
DOY.Array[,,i] <- terra::as.array(terra::rast(annualcrops[i], lyrs=DOY))
}
}
if(!is.na(PR)){
if(verbose){print(paste0("Creating Pixel Reliability array... ", Sys.time()))}
PR.Array <- array(dim=c(dim(terra::as.array(terra::rast(annualcrops[1])))[1],
dim(terra::as.array(terra::rast(annualcrops[1])))[2],
length(annualcrops)),
data=NA)
for(i in 1:length(annualcrops)){
PR.Array[,,i] <- terra::as.array(terra::rast(annualcrops[i], lyrs=PR))
}
}
# Clean NDVI data based on VIQ and PR data
PRcleaner <- function(NDVISet, PRSet){
Data <- NDVISet
Meta <- PRSet
for (i in 1:length(Data)){
if (!is.na(Meta[i])){
if (Meta[i] > 1){
Data[i] <- NA
}
if(Meta[i] == -1){
Data[i] <- NA
}
}
}
return(Data)
}
VIQcleaner <- function(NDVISet, VIQSet){
Data <- NDVISet
Meta <- VIQSet
MODLAND.QA <- as.numeric(substring(Meta, 1, 2))
VI.Usefulness <- as.numeric(substring(Meta, 3, 6))
Aerosol.Quantity <- as.numeric(substring(Meta, 7, 8))
Adjacent.Clouds <- as.numeric(substring(Meta, 9))
BRDF.Correction <- as.numeric(substring(Meta, 10))
Mixed.Clouds <- as.numeric(substring(Meta, 11))
Land.Water.Flag <- as.numeric(substring(Meta, 12, 14))
Possible.Snow.Ice <- as.numeric(substring(Meta, 15))
Possible.Shadow <- as.numeric(substring(Meta, 16))
for (i in 1:length(Data)){
if (MODLAND.QA[i] > 10){
Data[i] <- NA
}
if (VI.Usefulness[i] > 1100){
Data[i] <- NA
}
}
return(Data)
}
if(verbose){print(paste0("Cleaning NDVI data... ", Sys.time()))}
#FinalNDVI.PR <- PRcleaner(NDVISet = NDVIcleaned, PRSet = PR.Array)
FinalNDVI.VIQ <- VIQcleaner(NDVISet = NDVIcleaned, VIQSet = VIQbinary)
FinalNDVI <- FinalNDVI.VIQ
## Prep your NDVI values for Phenex processing
## Begin by creating a 3-D array in which
## NDVI values will are sorted by date of composite creation
if(verbose){
print(paste0("Sorting NDVI data into complete annual space-time cube... ",
Sys.time()))
}
year_array <- array(data = NA, dim = c(nrow(FinalNDVI), ncol(FinalNDVI), 366))
for (p in 1:dim(FinalNDVI)[3]){
zindex <- p
for (i in 1:nrow(FinalNDVI)){
rowindex <- i
for (j in 1:ncol(FinalNDVI)){
colindex <- j
if (!is.na(DOY.Array[rowindex, colindex, zindex])){
if(DOY.Array[rowindex, colindex, zindex] > 0){
year_array[
rowindex, colindex, DOY.Array[
rowindex, colindex, zindex]] <-
FinalNDVI[rowindex,colindex, zindex]
}
}
}
}
}
if(verbose){
print(paste0("Converting annual space-time cube into phenoscape... ", Sys.time()))
print("NDVI.Array pixel example")
print(NDVI.Array[12,15,])
print("VIQ Processed NDVI pixel example")
print(FinalNDVI.VIQ[12,15,])
print("Composite Day of Year pixel example")
print(DOY.Array[12,15,])
print("Annual space-time cube pixel example")
print(year_array[12,15,])
}
## Create a matrix to be filled with phenology data
phenodump <- array(data=NA, dim = c(nrow(FinalNDVI), ncol(FinalNDVI)))
## We want to apply the phenex functions through
## time in the FinalNDVI array;
## Meaning use plyr::aaply in dimensions c(1,2)
NDVIFunction <- function(a){
pixelID <- a
ndvi <- phenex::modelNDVI((pixelID/10000),
correction = "none",
method = "DLogistic",
year.int=as.integer(year))[[1]]
springtime <- phenex::phenoPhase(ndvi,
phase = phase,
method = "local",
threshold = threshold)[[1]]
return(springtime)
}
phenodump <- plyr::aaply(.data = year_array,
.fun = NDVIFunction,
.margins = c(1,2),
.progress = "text")
if(verbose){
print("Phenoscape produced...")
print(paste0("Converting phenoscape to raster... ", Sys.time()))
print(class(phenodump))
print(dim(phenodump))
}
## Turn final array back into a raster based on original raster properties
templateRast <- terra::rast(annualcrops[[1]],lyrs=1)
phenoscape <- terra::rast(x = phenodump)
terra::crs(phenoscape) <- terra::crs(templateRast)
terra::ext(phenoscape) <- terra::ext(templateRast)
# plot(phenoscape)
}
if(PhenoFactor == "Snow"){
if(verbose){print("Snowmelt option selected...")}
if(is.null(SnowExtent)){
stop("Crucial Dataset missing (have you included a SnowExtent layer?)")
}
Extent.Stack <- terra::rast(annualcrops,
lyrs = SnowExtent)
Extent.Array <- terra::as.array(Extent.Stack)
Complete.Extent.Array <- array(dim = c(dim(Extent.Array)[1],
dim(Extent.Array)[2],
366))
#Sort Extent.Array into daily slices (rather than 8-day slices)
for(p in (1:dim(Extent.Array)[3])){
zindex <- p
zfinal <- p*8-7
for(q in (1:nrow(Extent.Array))){
rowindex <- q
for(r in (1:ncol(Extent.Array))){
colindex <- r
Complete.Extent.Array[rowindex,colindex,zfinal] <-
Extent.Array[rowindex,colindex,zindex]
}
}
}
# Create a phenodump array to deposit snowmelt phenology data
phenodump.extent <- array(data=NA,
dim = c(nrow(Complete.Extent.Array),
ncol(Complete.Extent.Array)))
# Fill in phenodump with the first date at which x == 25, that is,
# the first time slice when pixels are detected as completely snow-free
phenodump.extent <- plyr::aaply(.data = Complete.Extent.Array,
.margins = c(1,2),
.fun = function(x) min(which(x == 25)))
templateRast <- terra::rast(annualcrops[[1]],lyrs=1)
phenoscape <- terra::rast(x = phenodump.extent)
terra::crs(phenoscape) <- terra::crs(templateRast)
terra::ext(phenoscape) <- terra::ext(templateRast)
}
if(verbose) {print(Sys.time())}
return(phenoscape)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.