编程小白近一星期的研究，走了不少弯路，一开始想要不自己写一下如何计算SPEI，实在是太高看自己编程小白的能力了，接下来自然不了了之了，后来后来，还是直接将R包拿来用吧。。
上手之前还是来看一下spei计算原理，话不多说，直接copy，copy图片来自于“1961–2015年新疆区域SPEI干旱指数数据集”介绍，http://www.csdata.org/p/217/
image
image
看完原理之后就基本可以调用R包了，首先需要安装SPEI包，install.package("SPEI")就ok了，Rstidio自动进行下载并安装。R包的spei包非常简便，可以直接拿来调用，建议改过程还是要参看官方帮助文档，讲解的很清楚也有简单实例。
其中计算PET潜在蒸散时，R包提供三种计算方法，分别是thornthwaite(),hargreaves(),penman().其中第一种方法最为简单，仅需要温度和纬度（单位：°）作为输入数据即可。以此准备数据：
数据准备
（1）温度数据，每一个站点按列排列，依次为开始年份的12个月，后一年的12个月，以此类推，简单来说就是一个站点所有年份一列数据，类似这个样子，格式是要保存为以空格分割的csv文件或txt文件，下同
image
（2）纬度数据
image
（3）降水数据（格式同温度数据格式）
准备完毕数据之后，就开始直接计算啦，话不多说，直接上代码. 该代码不是直接输出最终结果spei值，有输出的每一步结果，需要键入的路径可能比较多，小白的代码可能比较笨拙（为了自己好理解），仅供参考。
library(SPEI)
# ===========读入数据
# 降水数据
prec = read.table('G:\\test\\test3\\prec.txt',header = T)
# 温度数据
temp = read.table('G:\\test\\test3\\temp.txt',header = T)
# 每一个站点对应的纬度数据
lat = read.table('G:\\test\\test3\\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[,i] = thornthwaite(temp[,i],lat[,i])
PET_1 = cbind(PET_1,PET[,i])
}
write.table(PET_1,"G:\\test\\test3\\PET\\PET_station1.csv",sep = ',' ,row.names = FALSE)
# 到此保存的csv的列名均相同，需进行合并后手动更改
# ==========================================计算水平衡
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的列名均相同，需进行合并后手动更改
# ========================================= 计算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的输出格式没有明白，下列保存方法出错
# 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 = "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]，需进行合并后手动更改
}
# # 将采用上述方法得到的独立的站点spei+编号的csv文件按列合并，编号要手动更改，推荐用total commander，
# 比如上述得到的文件夹命名序号为1,2.....800,程序读进的顺序则会是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)
网友评论