怎么计算整个区域的SPEI啊

cahlogo

library(SPEI)
library(RNetCDF)
#read data
nc1 <- open.nc("C:\\Users\\cahlo\\OneDrive\\Desktop\\CRUTS\\cru_ts4.02degree1_tmp.nc")
nc2 <- open.nc("C:\\Users\\cahlo\\OneDrive\\Desktop\\CRUTS\\cru_ts4.02degree1_pre.nc")

tmp <-var.get.nc(nc1,"tmp")
precip <-var.get.nc(nc2,"pre")
lat <- var.get.nc( nc1, "lat")
lon <- var.get.nc( nc1, "lon")
time <-var.get.nc( nc1, "time")
spei3 <- array(NaN,c(360,180,length(time)-12))

tmp[is.na(tmp)]<-9.9692e+36
precip[is.na(precip)]<-9.9692e+36
for(i in 1:length(lon)&j in 1:length(lat)){

  BAL<- precip[i,j,c(1,1404)]-thornthwaite(tmp[i,j,c(1,1404)],lat=lat[j])

}
spei3<- spei(BAL, 3)   


出现错误Error in embed(acu, scale) : wrong embedding dimension
这个程序实在不会写，读取的是nc数据，求大神指点（之前再MATLAB计算了SPI，但是在没找到matlab的SPEI function，只能转战R。

dongyi

请问您是怎么用matlab计算整个区域的SPI呢

yulan

请问楼主做出来了吗

杰西卡·芬

楼主请问做出来了没有，求代码。我也是R语言做不出来。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)

# pdsi index
# mainland China. (GB/T 20481-2017)
options(PDSI.coe.K1.1 = 1.6)
options(PDSI.coe.K1.3 = 0.4)
options(PDSI.coe.K2 = 16.84)
options(PDSI.p = 0.755)
options(PDSI.q = 1/1.63)

for (i in 1:nlat) {
  for (j in 1:nlon) {
    ds[, "lat"] <- lat[i]
    ds[, "lon"] <- lon[j]
    # Computation of potential evapotranspiration
    ds[, "SPEI"] <- tryCatch({
      as.numberic(spei(BAL[j,i,], 1))
    }, 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-NorESM-1400-2000.nc", vars = indices, force_v4 = TRUE)
ncvar_put(ncnew, varid = "SPEI", vals = dd$SPEI)
nc_close(ncnew)