[R-sig-Geo] How can I apply function using stackApply or Calc, and save the results as raster object in R?

[Isaque Daniel]

Hi everyone, 
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"). 




My code dont save the results in two raster layers. 





Attached are the data (one small portion) 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




------------------------------------------------------------------------------------------------------------------Agronomist engineerMaster in Remote Sensing - National  Institute for Space Research (INPE) - BrazilPHD Student in Transport - Brasília University (UNB)  		 	   		  
	[[alternative HTML version deleted]]