格点数据算spei出错。请大神求助。

杰西卡·芬

有大神做过格点的spei吗？求代码。R官网只有一条序列的例子，我改成格点的但失败了。（其中能算出来pet和BAL），但是最后一步用spei函数算出来的都是缺测NA。我的代码如下：
setwd("D:/T_model/pet")
rm(list = ls())
library(lmomco)
library(SPEI)
library(scPDSI)
library(ncdf4)
library(plyr)
library(data.table)
library(tcltk)
diri1 <- "D:/T_model/pet/apt1400-2000.nc"
diri2 <- "D:/T_model/pet/cnz1400-2000.nc"
diri3 <- "D:/T_model/pet/cntsun1400-2000.nc"
diri4 <- "D:/T_model/pet/cnu1400-2000.nc"
diri5 <- "D:/T_model/pet/PM_PET-NorESM-1400-2000.nc"
diri6 <- "D:/T_model/pet/PM_BAL-NorESM-1400-2000.nc"
#diri1 <- "D:/T_model/pet/pr_1400-2000.nc"
#diri2 <- "D:/T_model/pet/t_1400-2000.nc"
#diri3 <- "D:/T_model/pet/tsun.nc"
#diri4 <- "D:/T_model/pet/z.nc"
#diri5 <- "D:/T_model/pet/AWND1000_1400-2000.nc"
##diri2 <- "D:/re-cmip5/past1000/FGOALS-gl/pr_Amon_FGOALS-gl_past1000_r1i1p1_160001-179912.nc"
#f1 <- nc_open(diri1)
#f2 <- nc_open(diri2)
#f3 <- nc_open(diri3)
#f4 <- nc_open(diri4)
#f5 <- nc_open(diri5)
##f2 <- nc_open(diri2)
f1 <- nc_open(diri1)
f2 <- nc_open(diri2)
f3 <- nc_open(diri3)
f4 <- nc_open(diri4)
f5 <- nc_open(diri5)
f6 <- nc_open(diri6)
#tmin <- ncvar_get(f2, varid = "tmin")
#tmin1 <- tmin-273.15
#tmax <- ncvar_get(f2, varid = "tmax")
#tmax1 <- tmax-273.15
#Tsun <- ncvar_get(f3, varid = "tsun")
#tsun1 <- Tsun*30.41
#z <- ncvar_get(f4, varid = "z")
#AWND <- ncvar_get(f5, varid = "AWND")
#pr <- ncvar_get(f1, varid = "pr")
#pp <- pr*24*30*3600
#pr1 <- pp*1000
tmin <- ncvar_get(f1, varid = "tmin")
tmin1 <- tmin-273.15
tmax <- ncvar_get(f1, varid = "tmax")
tmax1 <- tmax-273.15
z <- ncvar_get(f2, varid = "z")
Tsun <- ncvar_get(f3, varid = "tsun")
tsun1 <- Tsun*30.41
u2 <- ncvar_get(f4, varid = "u2")
AWND <- u2*3.6
pet <- ncvar_get(f5, varid = "pet")
BAL <- ncvar_get(f6, varid = "BAL")
pr <- ncvar_get(f1, varid = "pr")
pp <- pr*24*30*3600
pr1 <- pp*1000
##get the dimension
###close file and useless variables
lon <- ncvar_get(f1, varid = "lon")
lat <- ncvar_get(f1, varid = "lat")
##nc_close(f2)
rm(tmin)
rm(tmax)
rm(Tsun)
rm(u2)
rm(pr)
rm(pp)

###define time and empty array
time <- seq(from = as.Date("1400-02-01"), to = as.Date("2001-01-01"), by = "month")
time <- strsplit(as.character(time), split = "-")
time <- as.numeric(paste0(sapply(time, "[", 1), sapply(time, "[", 2)))
nlon <- length(lon)
nlat <- length(lat)
ntime <- length(time)
dd <- data.frame(time = numeric(), lat = numeric(), lon = numeric(), SPEI = numeric())
ds <- data.frame(time = time,
                 lat  = rep(NA, ntime),
                 lon  = rep(NA, ntime),
                 SPEI = rep(NA, ntime))
pb <- tkProgressBar(title = "cal-PDSI 1.0", label = "It has completed  %", min = 0, max = 100)

for (i in 1:nlat) {
  for (j in 1:nlon) {
    ds[, "lat"] <- lat[i]
    ds[, "lon"] <- lon[j]

    ds[, "SPEI"] <- tryCatch({
      as.numberic(spei(BAL[1,i,], 12))
    }, error = function(e) return(NA))
    dd <- rbind(dd, ds) # merge file
    info <- sprintf("It has completed %d%%", round(i*100/nlat))
    setTkProgressBar(pb, value = i*100/nlat, label = info)
  }
}
close(pb)
dd <- plyr::arrange(dd, time, lat, lon)
dd[is.na(dd$SPEI), "SPEI"] <- -999.9

lon <- ncvar_get(f1, varid = "lon")
lat <- ncvar_get(f1, varid = "lat")
#time <- ncvar_get(f1, varid = "time")
nc_close(f1)
lon <- ncdim_def(name = "lon",   units = "degreesE", vals = lon, longname = "Longitude_East")
lat <- ncdim_def(name = "lat",   units = "degreesN", vals = lat, longname = "Latitude_North")
time <- ncdim_def(name = "time", units = "days since 1001-01-01",
                  vals = time, unlim = TRUE, calendar = "noleap")
# define variables
indices <- ncvar_def(name = "SPEI", units = "", dim = list(lon, lat, time), missval = -999.9)
# open/create .nc file

ncnew <- nc_create(filename = "PM_SPEI-1400-2000.nc", vars = indices, force_v4 = TRUE)
ncvar_put(ncnew, varid = "SPEI", vals = dd$SPEI)
nc_close(ncnew)