I have been working in a list of function to process MODIS vegetation index timeseries and return crop cycle paramenters. It's for PhD thesis and I'll process all South America, so I need to use all forms to reduce the time to process.
I found the rasterEngine from the spatial.tools package to use paralell processing to speed up the process. but, before that I prepare some of the functions to compute by pixel over the raster stack my variables measures.
I develop the functions which will produce 7 different outputs, I tried tu use my function "CropAnalysis" to compute by each pixel, in the code in the post I try to save a raster brick with 2 layers (each one with one of variables produced by the function "CropAnalysis").
I edited the code, but couldn't solve the problem when I run the process.
Attached are the data (one small portion of data) and the code, any idea?
My data: Modis stack https://www.dropbox.com/s/uesgzv125e3v3e6/stackimagesNDVI.tif?dl=0
My code:
library(stringr)
library(rgdal)
library(raster)
# loading the data
limit <- 3000 # minimum value betweem maximum and minimum to be crop
ndates <- 2 # time difference between maximum and minimum to be crop
min_diff <- 3000 # threshold for the maximum value (this is the minimum value to test)
min_val <- 1500 # minimum value for the minimum pixel values be trustfull (threshold)
max_phase_duration <- 7 # the maximum interval over the maximum value between the two adjacent minimum values
number_of_crop_cycles <- 3 # definition of number of crop cycles per croo year
imgStacked <- brick('stackimagesNDVI.tif')
CropAnalysis <- function (pixel, ...){
pixel <- as.vector(pixel)
# test : if is No data the return is
if (is.na(pixel)) {-1}
else{
# delta (valor i - valor i+1)
delta <- pixel[2:length(pixel)] - pixel[1:(length(pixel)-1)]
# maximum and minimum point
ptma<-NULL
ptmi <- NULL
# verifing why the first time point is not signed???? T or F
if (pixel[2] > pixel [1]) {ptmi <- 1}
if (pixel[2] < pixel [1]) {ptma <- 1}
# computing the slope of the line change from positive to negative
for (j in 1:(length(delta)-1))
{
if (delta[j]>0 && delta[j+1]<0 )
{
ptma<- c(ptma,j+1) # point of maximum
}
if (delta[j]<0 && delta[j+1]>0)
{
ptmi<- c(ptmi,j+1) # point of minimum
}
}
# verifing why the first time point is not signed???? T or F
if (pixel[(length(pixel))] > pixel [(length(pixel)-1)]) {ptma <- c(ptma,length(pixel))}
if (pixel[(length(pixel))] < pixel [(length(pixel)-1)]) {ptmi <- c(ptmi,length(pixel))}
# variables for save the measures for crop cycle
max_points <- as.numeric(rep(0, number_of_crop_cycles)) # number of maximum peaks after test if is a crop pixel
length_max_period <- as.numeric(rep(NA, number_of_crop_cycles)) # variation of number of dates between the minimum points around of maximum point
max_valids <- NULL
# agricultural detection
for (j in 1:length(ptma))
{
index <- ptma[j]
# logical tests to verify the presence of crop
# from each maximum value, check if:
# 1st - the maximum position had the before minimum value far or equal than "ndatas" variable
# 2nd - the maximum position had the after minimum value far or equal than ndatas variable
# 3th - the value of maximum is equal or great than "val_min" variable (threshold)
# 4th - the difference between the maximum value and the two minimum values (in the "ndates") distance is equal or bigher than "limit" variable (threshold value of increase Vegetation index)
# 5th - the minimum values bigher tha minimum limit variable
# 6th - check to exclude sugarcane from anual crop cicle
if(!is.na(((ptmi[ptmi < index][length(ptmi[ptmi < index])]+ndates) <= index && # 1st test
index <= (ptmi[ptmi > index][1]-ndates)) && # 2sd test
(pixel[index] >= limit) && # 3th test
((pixel[index]-pixel[ptmi[ptmi < index]][length(pixel[ptmi[ptmi < index]])] >= min_diff) && (pixel[index]-pixel[ptmi[ptmi > index]][1] > min_diff)) && # 4th test
(pixel[ptmi[ptmi < index]][length(pixel[ptmi[ptmi < index]])] && pixel[ptmi[ptmi > index]][1] >= min_val) && # 5th
((ptmi[ptmi < index][length(ptmi[ptmi < index])] <= index-(max_phase_duration-3) && index-(max_phase_duration-3)>= 1) | (ptmi[ptmi > index][1] >= index+(max_phase_duration-3))))) # 6th
{
# computing the valid maximum values to avoid the "fake" crop pattern (small difference between min and max) and using this "position_data" to save the values in the vectors in the right order
max_valids <- c(max_valids, index)
position_data <- which(max_valids==index)
# saving the points of maximum per pixel over the time series
max_points[position_data] <- index
# calculating the crop cycle length
length_max_period[position_data] <- (index-ptmi[ptmi < index][length(ptmi[ptmi < index])])+(ptmi[ptmi > index][1]-index)
}
}
# replacing the NA data (NA is the default value and show possible cropseasons whitout crops)
#max_points[is.na(max_points)]<-0
# join the values in a unique number: i.e = c(5,16, 0) -> 99051600 ( 99 = to avoid the difference of length of pixel value in cases of numbers lower than 10; all valid number using flag 0)
max_points <- as.integer(paste('99',paste(formatC(max_points, flag=0, digits = 1,format = 'd'),collapse = ''),sep=""))
length_max_period <- as.integer(paste('99',paste(formatC(length_max_period, flag=0, digits = 1,format = 'd'),collapse = '')),sep="")
}
}
. # using stackApply
data_process <- stackApply(imgStacked, indices=c(rep(1,nlayers(imgStacked)),rep(2,nlayers(imgStacked))), fun=CropAnalysis)
The error message:
Error in length_max_period[position_data] <- (index - ptmi[ptmi < index [length(ptmi[ptmi < : replacement has length zero
In addition: Warning messages:
1: In stackApply(imgStacked, indices = c(rep(1, nlayers(imgStacked)), : number of items to replace is not a multiple of replacement length
2: In if (is.na(pixel)) { : the condition has length > 1 and only the first element will be used
# using calc
data_process<-calc(x=imgStacked, fun=CropAnalysis, forcefun=TRUE, forceapply=TRUE)
The error message:
Error in colnames<-(*tmp*, value = "layer") : length of 'dimnames' [2] not equal to array extent
In addition: Warning messages:
1: In if (is.na(pixel)) { : the condition has length > 1 and only the first element will be used
2: In if (is.na(pixel)) { :the condition has length > 1 and only the first element will be used
3: In fun(tstdat) : NAs introduced by coercion to integer range
4: In fun(tstdat) : NAs introduced by coercion
5: In if (is.na(pixel)) { :the condition has length > 1 and only the first element will be used
6: In fun(x) : NAs introduced by coercion to integer range
7: In fun(x) : NAs introduced by coercion
8: In matrix(values, nrow = ncell(x), ncol = nlayers(x)) : data length exceeds size of matrix
