R计算SPEI指数

library(SPEI)

# ===========读入数据

# 降水数据
prec = read.table('D:\\station\\prec.txt',header = T) 
# 温度数据
temp = read.table('D:\\station\\temp.txt',header = T)
# 每一个站点对应的纬度数据
lat = read.table('D:\\station\\lat.txt',header = T)
# 输出路径

print(ncol(temp))
print(ncol(prec))

# ==========================================计算PET 本代码采用thornthwaite方法计算，只需要温度数据
# 利用cbind函数将每一列连接，就需要初始化一列
PET_1 = thornthwaite(temp[,1],lat[,1])
for (i in 2:ncol(temp))
{
  PET = thornthwaite(temp[,i],lat[,i])
  PET_1 = cbind(PET_1,PET)
}
write.table(PET_1,"D:\\station\\PET_station1.csv",sep = ',' ,row.names = FALSE)
# 到此保存的csv的列名均相同，需进行合并后手动更改
# ==========================================计算水平衡
PET = read.table("D:\\station\\PET_station.txt",header = T)
balance = prec[,1]-PET[,1]
for (j in 2:ncol(prec))
{
  balance1 = prec[,j]-PET[,j]
  balance = cbind(balance,balance1)
}
write.table(balance,"D:\\station\\waterbalance_station.csv",sep = ",",row.names = FALSE) 
# 到此保存的csv的列名均相同，需进行合并后手动更改
# ========================================= 计算SPEI
library(SPEI)
water_bal = read.table("D:\\station\\waterbalance_station.txt",header = T)
spei12 = spei(water_bal[,1], 6)
print(spei1)
# write.csv(spei12$fitted,"D:\\station\\spei1.csv",row.names = F)
print(ncol(water_bal))
# # 因为spei的输出格式没有明白，下列保存方法出错
# for(k in 2:ncol(water_bal))
# {
#   spei1 = spei(water_bal[,k], 12)
#   spei12 = cbind(spei12,spei1)
#   write.table(spei12$fitted,"D:\\spei1.csv",row.names = F)
# }

# # 采用该方法将每一个站点的spei单独输出成独立的csv文件，较麻烦，后续需要将独立的csv文件合并
output = "D:\\station\\spei"
for(k in 1:ncol(water_bal))
{
  spei1 = spei(water_bal[,k], 6)
  output1 = paste(k,".csv",sep='')
  write.table(spei1$fitted,paste(output,output1,sep=''),row.names = F)  # 到此所有csv的列名均为water_bal[,k]，需进行合并后手动更改
}

# # 将采用上述方法得到的独立的站点spei+编号的csv文件按列合并，编号要手动更改，推荐用total commander，
#   比如上述得到的文件夹命名序号为1,2.....800,程序读进的顺序则会是1,10,11,12,13...19,100,101这样的顺序，我将名称均改为三位的，如001，002,003
input_dir = "D:\\station\\spi\\"
file = list.files(input_dir)
dir = paste(input_dir,file,sep="")
n = length(dir)
merge_data = read.csv(dir[1],header = T)
for (i in 2:n){
  new.data = read.csv(dir[i], header=T)
  merge_data = cbind(merge_data,new.data)
}
write.csv(merge_data,file = "D:\\spei1.csv",row.names=F)