[R-sig-Geo] How can I apply function using stackApply or Calc, and save the results as raster object in R?
Isaque Daniel
isaquedanielre at hotmail.com
Thu Apr 28 18:34:48 CEST 2016
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]]
More information about the R-sig-Geo
mailing list