編程小白近一星期的研究既绩，走了不少彎路，一開始想要不自己寫一下如何計算SPEI私杜，實在是太高看自己編程小白的能力了救欧，接下來自然不了了之了笆怠，后來后來，還是直接將R包拿來用吧蹬刷。频丘。
上手之前還是來看一下spei計算原理，話不多說迂卢，直接copy桐汤，copy圖片來自于“1961–2015年新疆區(qū)域SPEI干旱指數(shù)數(shù)據(jù)集”介紹怔毛，http://www.csdata.org/p/217/
看完原理之后就基本可以調(diào)用R包了馆截，首先需要安裝SPEI包蜡娶，install.package("SPEI")就ok了，Rstidio自動進(jìn)行下載并安裝幕随。R包的spei包非常簡便宿接，可以直接拿來調(diào)用睦霎，建議改過程還是要參看官方幫助文檔，講解的很清楚也有簡單實例蛤高。
其中計算PET潛在蒸散時，R包提供三種計算方法塞绿，分別是thornthwaite(),hargreaves(),penman().其中第一種方法最為簡單异吻，僅需要溫度和緯度（單位：°）作為輸入數(shù)據(jù)即可喜庞。以此準(zhǔn)備數(shù)據(jù)：
數(shù)據(jù)準(zhǔn)備
（1）溫度數(shù)據(jù)延都，每一個站點按列排列，依次為開始年份的12個月春宣，后一年的12個月嫉你，以此類推，簡單來說就是一個站點所有年份一列數(shù)據(jù)嚷辅，類似這個樣子簸搞，格式是要保存為以空格分割的csv文件或txt文件，下同
（2）緯度數(shù)據(jù)
（3）降水?dāng)?shù)據(jù)（格式同溫度數(shù)據(jù)格式）
準(zhǔn)備完畢數(shù)據(jù)之后，就開始直接計算啦刑然，話不多說泼掠，直接上代碼. 該代碼不是直接輸出最終結(jié)果spei值，有輸出的每一步結(jié)果择镇，需要鍵入的路徑可能比較多腻豌，小白的代碼可能比較笨拙（為了自己好理解）叹谁，僅供參考。
library(SPEI)
# ===========讀入數(shù)據(jù)
# 降水?dāng)?shù)據(jù)
prec = read.table('G:\\test\\test3\\prec.txt',header = T)
# 溫度數(shù)據(jù)
temp = read.table('G:\\test\\test3\\temp.txt',header = T)
# 每一個站點對應(yīng)的緯度數(shù)據(jù)
lat = read.table('G:\\test\\test3\\lat.txt',header = T)
# 輸出路徑
print(ncol(temp))
print(ncol(prec))
# ==========================================計算PET 本代碼采用thornthwaite方法計算订框，只需要溫度數(shù)據(jù)
# 利用cbind函數(shù)將每一列連接兜叨，就需要初始化一列
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,"G:\\test\\test3\\PET\\PET_station1.csv",sep = ',' ,row.names = FALSE)
# 到此保存的csv的列名均相同国旷，需進(jìn)行合并后手動更改
# ==========================================計算水平衡
PET = read.table("G:\\test\\test3\\PET\\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,"G:\\test\\test3\\PET\\waterbalance_station.csv",sep = ",",row.names = FALSE)
# 到此保存的csv的列名均相同履羞，需進(jìn)行合并后手動更改
# ========================================= 計算SPEI
library(SPEI)
water_bal = read.table("G:\\test\\test3\\PET\\waterbalance_station.txt",header = T)
spei12 = spei(water_bal[,1], 12)
print(spei12)
# write.csv(spei12$fitted,"D:\\spei1.csv",row.names = F)
print(ncol(water_bal))
# # 因為spei的輸出格式?jīng)]有明白屡久，下列保存方法出錯
# 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文件被环，較麻煩，后續(xù)需要將獨立的csv文件合并
output = "G:\\test\\test3\\PET\\spei\\spei"
for(k in 1:ncol(water_bal))
{
spei1 = spei(water_bal[,k], 12)
output1 = paste(k,".csv",sep='')
write.table(spei1$fitted,paste(output,output1,sep=''),row.names = F) # 到此所有csv的列名均為water_bal[,k]浸锨，需進(jìn)行合并后手動更改
}
# # 將采用上述方法得到的獨立的站點spei+編號的csv文件按列合并柱搜，編號要手動更改剥险，推薦用total commander冯凹，
# 比如上述得到的文件夾命名序號為1,2.....800,程序讀進(jìn)的順序則會是1,10,11,12,13...19,100,101這樣的順序，我將名稱均改為三位的炒嘲，如001宇姚，002,003
input_dir = "G:\\test\\test3\\PET\\spei\\"
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:\\spei.csv",row.names=F)
