MODIS 植生指数の時系列を処理し、作物サイクルのパラメーターを返す関数のリストで作業してきました。これは博士論文用で、南米全体を処理するので、処理時間を短縮するためにすべてのフォームを使用する必要があります。
並列処理を使用してプロセスを高速化するために、spatial.tools パッケージの rasterEngine を見つけました。ただし、その前に、変数が測定するラスター スタックをピクセル単位で計算する関数をいくつか準備します。
私は7つの異なる出力を生成する関数を開発し、関数「CropAnalysis」を使用して各ピクセルごとに計算しようとしました。投稿のコードで、2つのレイヤーでラスターブリックを保存しようとしました(それぞれに生成された変数の1つがあります)関数「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 < index [length(ptmi[ptmi < : 置換の長さはゼロです)
さらに: 警告メッセージ:
1: stackApply(imgStacked, インデックス = c(rep(1, nlayers(imgStacked))) 内: 置換するアイテムの数が置換長の倍数ではありません
2: in if (is.na(pixel)) { : 条件の長さ > 1 で、最初の要素のみが使用されます
# using calc
data_process<-calc(x=imgStacked, fun=CropAnalysis, forcefun=TRUE, forceapply=TRUE)
エラーメッセージ:
colnames<-( 、値 = "layer") のエラー*tmp*: 'dimnames' [2] の長さが配列範囲と等しくありません
さらに: 警告メッセージ:
1: in if (is.na(pixel)) { : 条件の長さ > 1 で、最初の要素のみが使用されます
2: if (is.na(pixel)) { : 条件の長さ > 1 で、最初の要素のみが使用されます
3: in fun(tstdat) : 整数範囲への強制によって導入された NA
4: In fun(tstdat) : 強制によって導入された NA
5: if (is.na(pixel)) { : 条件の長さ > 1 で、最初の要素のみが使用されます
6: in fun(x) : 整数範囲への強制によって導入された NA
7: In fun(x) : 強制によって導入された NA
8: matrix(values, nrow = ncell(x), ncol = nlayers(x)) 内: データ長が行列のサイズを超えています