0

我一直在处理 MODIS 植被指数时间序列和返回作物周期参数的函数列表。这是博士论文，我会处理整个南美洲，所以我需要使用所有表格来减少处理时间。

我从 spatial.tools 包中找到了 rasterEngine 来使用并行处理来加速这个过程。但是，在此之前，我准备了一些函数，以便在我的变量测量的栅格堆栈上按像素计算。

我开发了将产生 7 个不同输出的函数，我尝试使用我的函数“CropAnalysis”来计算每个像素，在帖子中的代码中我尝试保存一个有 2 层的光栅砖（每一个都有一个产生的变量通过功能“作物分析”）。

我编辑了代码，但在运行该过程时无法解决问题。

附上数据（一小部分数据）和代码，知道吗？

我的数据：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, indices = c(rep(1, nlayers(imgStacked)), : 要替换的项目数不是替换长度的倍数

2: 在 if (is.na(pixel)) { : 条件长度 > 1 并且只使用第一个元素

    # using calc
    data_process<-calc(x=imgStacked, fun=CropAnalysis, forcefun=TRUE, forceapply=TRUE)

错误信息：

colnames<-( , value = "layer") 中的错误*tmp*：'dimnames' [2] 的长度不等于数组范围

另外：警告信息：

1: 在 if (is.na(pixel)) { : 条件长度 > 1 并且只使用第一个元素

2: 在 if (is.na(pixel)) { : 条件长度 > 1 并且只使用第一个元素

3：在 fun(tstdat) 中：强制引入整数范围的 NA

4：在 fun(tstdat) 中：强制引入的 NA

5: 在 if (is.na(pixel)) { : 条件长度 > 1 并且只使用第一个元素

6：在 fun(x) 中：通过强制引入整数范围的 NA

7：在 fun(x) 中：强制引入的 NA

8：在矩阵中（值，nrow = ncell（x），ncol = nlayers（x））：数据长度超过矩阵的大小

4

0 回答 0