-
Notifications
You must be signed in to change notification settings - Fork 0
/
module1_rainfall.R
119 lines (93 loc) · 3.21 KB
/
module1_rainfall.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#mean(x)
#median(x)
#sd(x) #standard deviation
#hist(x,main="Histogram of observed data")
#z<-(x-mean(x))/sd(x) ##standardized data
#hist(z,main="Histogram of observed data")
#plot(density(x),main="Density estimate of data")
#plot(ecdf(x),main="Empirical cumulative distribution function")
#alfa<-(median(x)^2)/(sd(x)^2)
#beta<-(sd(x)^2)/(median(x)^2)
################################################
#par(mfrow = c(2, 2))
#denscomp(fg)
#qqcomp(fg)
#cdfcomp(fg)
#ppcomp(fg)
#xlab="Return period (years)",
#ylab=expression(paste("Discharge [",m^3,"/s]")),
#log="x")
##################################################
cum_3<-read.csv(file = "cumulative_3_months.csv")
Target<-"DJF"
x<-cum_3[,which(colnames(cum_3)==Target)] #extract all cumulative data in the target 3-month
#fit gamma dist. with MLE method:(fitdistrplus package)
library(fitdistrplus) ## loading package fitdistplus
#find the zeros
x.1<-data.frame(x)
num.zer<-colSums(x.1==0.0)
m <- unname(num.zer)
n <- length(x)
if (num.zer==0) {fg<-fitdist(x,"gamma",method ="mle" )$estimate #finding gamma parameters (PDF)
par(mfrow = c(1, 2))
plot(ecdf(x),col="blue",
main="Empirical and Gamama CDF",
xlab=c("3-month precipitation (mm)",Target) ,
ylab="Cumulative probability")
shape<-fg["shape"]
rate<-fg["rate"]
G = pgamma(x, shape, rate)
points(x,G,pch=18)
x.gamma <- seq(min(x),max(x),1)
G.gamma <- pgamma(x.gamma,shape,rate)
lines(x.gamma,G.gamma,lwd=2)
legend("bottomright",legend = c("Empirical CDF","Gamma CDF"),
col=c("blue", "black"),bty = "n",lty=0:1, pch =c(19,18), pt.cex = 1,
inset = c(-0.3, 0.02), cex=0.56)
#standardizing the data to reach SPI
spi<-qnorm(G)
plot(spi,G,pch=18,
main="Standard Normal CDF",
ylim=c(0,1), # y axis limits
xlab=c("SPI",Target) ,
ylab=""
)
G.gamma.nor <- qnorm(G.gamma)
lines(G.gamma.nor,G.gamma,lwd=2)
} else {
x_nonzero <- x
x_nonzero[x_nonzero == 0] <- NA #substituting zeros with NA
x_nonzero<-x_nonzero[!is.na(x_nonzero)] #removing Na from a vector
fg<-fitdist(x_nonzero,"gamma",method ="mle")$estimate #finding gamma parameters (PDF)
par(mfrow = c(1, 2))
plot(ecdf(x),col="blue",
main="Empirical and Gamama CDF",
xlab=c("3-month precipitation (mm)",Target) ,
ylab="Cumulative probability")
shape<-fg["shape"]
rate<-fg["rate"]
G = pgamma(x, shape, rate)
q <- m/n
H <- q+(1-q)*G
points(x,H,pch=18)
x.gamma <- seq(min(x),max(x),0.5)
G.gamma <- pgamma(x.gamma,shape,rate)
H.gamma <- q+(1-q)*G.gamma
lines(x.gamma,H.gamma,lwd=2)
legend("bottomright",legend = c("Empirical CDF","Gamma CDF"),
col=c("blue", "black"),bty = "n",lty=0:1, pch =c(19,18), pt.cex = 1,
inset = c(-0.3, 0.02), cex=0.56)
#standardizing the data to reach SPI
spi<-qnorm(H)
plot(spi,H,pch=18,
main="Standard Normal CDF",
ylim=c(0,1), # y axis limits
xlab=c("SPI",Target) ,
ylab=""
)
H.gamma.nor <- qnorm(H.gamma)
lines(H.gamma.nor,H.gamma,lwd=2)
}
shape;rate #gamma parameters
#library(xlsx)
#write.xlsx(spi,"DJF.xlsx")