Как я могу применить функцию с помощью stackApply или Calc в rasterEngine и сохранить результаты как растровый объект в R? # отредактировано
Я работал над списком функций для обработки временных рядов индекса вегетации MODIS и параметров возврата цикла урожая. Это для докторской диссертации, и я обработаю всю Южную Америку, поэтому мне нужно использовать все формы, чтобы сократить время на обработку.
Я нашел rasterEngine из пакета atial.tools, чтобы использовать обработку paralell для ускорения процесса. но перед этим я подготавливаю некоторые функции для вычисления по пикселам над стеком растров меры моих переменных.
Я разрабатываю функции, которые будут генерировать 7 различных выходных данных, я пытался использовать мою функцию "CropAnalysis" для вычисления по каждому пикселю, в коде поста я пытаюсь сохранить растровый кирпич с 2 слоями (каждый с одной из произведенных переменных функцией "CropAnalysis").
Я редактировал код, но не смог решить проблему при запуске процесса.
Приложены данные (одна небольшая часть данных) и код, есть идеи?
Мои данные: стек Modis https://www.dropbox.com/s/uesgzv125e3v3e6/stackimagesNDVI.tif?dl=0
Мой код:
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="")
}
}
, # используя stackApply
data_process <- stackApply(imgStacked, indices=c(rep(1,nlayers(imgStacked)),rep(2,nlayers(imgStacked))), fun=CropAnalysis)
Сообщение об ошибке:
Ошибка в length_max_period[position_data] <- (index - ptmi [ptmi Дополнительно: Предупреждающие сообщения: 1: в stackApply(imgStacked, indices = c(rep(1, nlayers(imgStacked)),: количество заменяемых элементов не кратно длине замены 2: В if (is.na(pixel)) {: условие имеет длину> 1, и будет использоваться только первый элемент Сообщение об ошибке: Ошибка в Дополнительно: Предупреждающие сообщения: 1: в if (is.na(pixel)) {: условие имеет длину> 1 и будет использоваться только первый элемент 2: В if (is.na(pixel)) {: условие имеет длину> 1, и будет использоваться только первый элемент 3: В шутку (tstdat): NA вводятся путем принуждения в целочисленный диапазон 4: В шутку (tstdat): НС вводятся по принуждению 5: In if (is.na(pixel)) {: условие имеет длину> 1, и будет использоваться только первый элемент 6: В шутке (х): NA введены путем принуждения к целочисленному диапазону 7: В шутку (х): НС, введенные по принуждению 8: В матрице (значения nrow = ncell(x), ncol = nlayers(x)): длина данных превышает размер матрицы # using calc
data_process<-calc(x=imgStacked, fun=CropAnalysis, forcefun=TRUE, forceapply=TRUE)
colnames<-(*tmp*, value = "layer"): длина 'dimnames' [2] не равна экстенту массива