[交流] R语言潜在蒸散发代码已有7人参与

####设置工作目录
setwd("C://Users//36242//Desktop")
####读取及数据预选择
firstdata<-read.csv("2000年气象数据.csv")
View(firstdata)
seconddata<-firstdata[,c(2:8,12:15,18,23)]
##选取计算的站点
thirddata<-subset(seconddata,Plot=="57445")
##将原数据月、日两列排序
newdata<-thirddata[order(thirddata$month,thirddata$day),]
##修改日列数据为每天
newdata$day<-c(1:sum(!is.na(newdata$Plot)))
####计算实际水汽压及饱和水汽压（参数用到最高温、最低温、以及平均相对湿度，温度数据气象数据共享网单位为0.1摄氏度，湿度单位为%）
#最高温水汽压
emax<-0.611*exp(17.27*(newdata$maxTEM/10)/((newdata$maxTEM/10)+237.3))
emax
##最低温水汽压
emin<-0.611*exp(17.27*(newdata$minTEM/10)/((newdata$minTEM/10)+237.3))
emin
#饱和水汽压
ea<-(emax+emin)/2
ea
#实际水汽压
ed<-(newdata$meanRHU)/((50/emin)+(50/emax))
ed
####饱和水汽压与温度之间关系曲线斜率
Rate<-(4098*ea)/(((newdata$meanTEM/10)+237.2)^2)
Rate
####干湿表常数(参数用到所在站点的高程、温度)
##高程H处气压P
P<-101.3*(((293-(0.0065*newdata$Elevation/10))/293)^5.26)
P
##或者u可以取值为2.501
u<-2.501-(0.002361*newdata$meanTEM/10)
u
##干湿常数
gamma<-0.00163*P/u
gamma
####计算2米高处的风速（平均风速、风速测量的高度，一般为气象观测高度要求大于10m,此处选取10m，ln在R中表达为log）
U<-4.87*newdata$meanWIN/10/(log((67.8*10-5.42)))
U
####土壤热通量（相对净辐射，可取值为0）
####太阳净辐射
#太阳的磁偏角（所用参数为当年的天)
Delta<-0.409*sin((2*pi*newdata$day/365)-1.39)
Delta
#日-地相对距离
dr<-1+(0.0338*cos(2*pi*newdata$day/365))
dr
#lat代表站点所在的纬度(3302即为33度22分；一般的纬度指的是角度，弧度=角度*pi/180)
lat1<-round(newdata$X/100,digits=0)
lat<-((((newdata$X/100)-lat1)*100/60)+lat1)*pi/180
lat
#日落时角度（反余弦函数表达式为acos）
Omega<-((acos(-tan(lat[1])*tan(Delta))))
Omega
#碧空太阳辐射
Ra<-37.6*dr*((Omega*sin(lat[1])*sin(Delta))+(cos(lat[1])*cos(Delta)*sin(Omega)))
Ra
#最大天文日照时数
N<-24*Omega/pi
N
#净长波辐射（Tkx表达为摄氏度+273.6）
Rl<-2.45*(10^-9)*(0.1+(0.9*newdata$SSD/10/N))*(0.34-(0.14*sqrt(ed)))*(((max(newdata$maxTEM/10)+273.6)^4)+((min(newdata$minTEM/10)+273.6)^4))
Rl
#净短波辐射（日照时数单位为0.1h）
Rn<-0.77*(0.25+(0.5*newdata$SSD/10/N))*Ra
Rn
sum(Rn)
sum(Rl)
#净辐射量
R<-Rn-Rl
R
View(R)
##########潜在蒸散发
ET<-((0.408*Rate*R)+(gamma*U*(ea-ed)*900/(273+(newdata$meanTEM/10))))/(Rate+(gamma*(1+(0.34*U))))
ET
sum(ET)

2楼: Originally posted by Lllitian at 2018-06-25 20:12:14
真的好棒，可以加微信共同学习吗？我也在研究invest

3楼: Originally posted by 勾蒙蒙 at 2018-06-26 00:02:52
可以，你微信多少
...