-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhw_rnorm_hist_90p.R
200 lines (154 loc) · 9.82 KB
/
hw_rnorm_hist_90p.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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
setwd("//Volumes/Extended/INMET/allstations/")
dfm <- read.csv("percentile_occ.csv", header = TRUE, sep = ",")
p1 <- rnorm(length(dfm$d1_tx90), mean(dfm$d1_tx90, na.rm = TRUE), sd(dfm$d1_tx90, na.rm = TRUE))
p2 <- rnorm(length(dfm$d2_tx90), mean(dfm$d2_tx90, na.rm = TRUE), sd(dfm$d2_tx90, na.rm = TRUE))
p3 <- rnorm(length(dfm$d1_tn90), mean(dfm$d1_tn90, na.rm = TRUE), sd(dfm$d1_tn90, na.rm = TRUE))
p4 <- rnorm(length(dfm$d2_tn90), mean(dfm$d2_tn90, na.rm = TRUE), sd(dfm$d2_tn90, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_All_sumtx90p.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=20, col=rgb(0.5,0,0,0.4), xlab="Days with Tx >90p", ylim=c(0,60), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p2, prob=FALSE, breaks=10, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
quartz()
png(filename = paste("rnorm_All_sumtn90p.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=20, col=rgb(0.5,0,0,0.4), xlab="Days with Tn >90p", ylim=c(0,60), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p4, prob=FALSE, breaks=20, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
CW <- subset(dfm, Region == "CW") # D1 as 2007-2016 including both years - 10y
N <- subset(dfm, Region == "N") # D1 as 2007-2016 including both years - 10y
NE <- subset(dfm, Region == "NE") # D1 as 2007-2016 including both years - 10y
S <- subset(dfm, Region == "S") # D1 as 2007-2016 including both years - 10y
SE <- subset(dfm, Region == "SE") # D1 as 2007-2016 including both years - 10y
#--------------- Central western SKIPPED only 5 stations -----------------
p1 <- rnorm(length(NE$d1_tx90), mean(NE$d1_tx90, na.rm = TRUE), sd(NE$d1_tx90, na.rm = TRUE))
p2 <- rnorm(length(NE$d2_tx90), mean(NE$d2_tx90, na.rm = TRUE), sd(NE$d2_tx90, na.rm = TRUE))
p3 <- rnorm(length(NE$d1_tn90), mean(NE$d1_tn90, na.rm = TRUE), sd(NE$d1_tn90, na.rm = TRUE))
p4 <- rnorm(length(NE$d2_tn90), mean(NE$d2_tn90, na.rm = TRUE), sd(NE$d2_tn90, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_NE_sumtx90p.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tx >90p", ylim=c(0,25), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p2, prob=FALSE, breaks=10, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1,25, "Northeast region", cex=0.9, col="grey20", adj = 0)
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
quartz()
png(filename = paste("rnorm_NE_sumtn90p.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=20, col=rgb(0.5,0,0,0.4), xlab="Days with Tn >90p", ylim=c(0,25), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p4, prob=FALSE, breaks=20, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1,25, "Northeast region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- North -----------------
p1 <- rnorm(length(N$d1_tx90), mean(N$d1_tx90, na.rm = TRUE), sd(N$d1_tx90, na.rm = TRUE))
p2 <- rnorm(length(N$d2_tx90), mean(N$d2_tx90, na.rm = TRUE), sd(N$d2_tx90, na.rm = TRUE))
p3 <- rnorm(length(N$d1_tn90), mean(N$d1_tn90, na.rm = TRUE), sd(N$d1_tn90, na.rm = TRUE))
p4 <- rnorm(length(N$d2_tn90), mean(N$d2_tn90, na.rm = TRUE), sd(N$d2_tn90, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_N_sumtx90p.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tx >90p", ylim=c(0,10), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p2, prob=FALSE, breaks=10, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1,10, "North region", cex=0.9, col="grey20", adj = 0)
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
quartz()
png(filename = paste("rnorm_N_sumtn90p.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tn >90p", ylim=c(0,10), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p4, prob=FALSE, breaks=10, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1,10, "North region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- Southeast -----------------
p1 <- rnorm(length(SE$d1_tx90), mean(SE$d1_tx90, na.rm = TRUE), sd(SE$d1_tx90, na.rm = TRUE))
p2 <- rnorm(length(SE$d2_tx90), mean(SE$d2_tx90, na.rm = TRUE), sd(SE$d2_tx90, na.rm = TRUE))
p3 <- rnorm(length(SE$d1_tn90), mean(SE$d1_tn90, na.rm = TRUE), sd(SE$d1_tn90, na.rm = TRUE))
p4 <- rnorm(length(SE$d2_tn90), mean(SE$d2_tn90, na.rm = TRUE), sd(SE$d2_tn90, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_SE_sumtx90p.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tx >90p", ylim=c(0,30), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p2, prob=FALSE, breaks=10, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1,30, "Southeast region", cex=0.9, col="grey20", adj = 0)
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
quartz()
png(filename = paste("rnorm_SE_sumtn90p.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tn >90p", ylim=c(0,30), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p4, prob=FALSE, breaks=10, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1,30, "Southeast region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- South -----------------
p1 <- rnorm(length(S$d1_tx90), mean(S$d1_tx90, na.rm = TRUE), sd(S$d1_tx90, na.rm = TRUE))
p2 <- rnorm(length(S$d2_tx90), mean(S$d2_tx90, na.rm = TRUE), sd(S$d2_tx90, na.rm = TRUE))
p3 <- rnorm(length(S$d1_tn90), mean(S$d1_tn90, na.rm = TRUE), sd(S$d1_tn90, na.rm = TRUE))
p4 <- rnorm(length(S$d2_tn90), mean(S$d2_tn90, na.rm = TRUE), sd(S$d2_tn90, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_S_sumtx90p.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tx >90p", ylim=c(0,15), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p2, prob=FALSE, breaks=5, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1,15, "South region", cex=0.9, col="grey20", adj = 0)
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
quartz()
png(filename = paste("rnorm_S_sumtn90p.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Days with Tn >90p", ylim=c(0,15), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p4, prob=FALSE, breaks=5, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1,15, "South region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- Centralwestern -----------------
p1 <- rnorm(length(CW$d1_tx90), mean(CW$d1_tx90, na.rm = TRUE), sd(CW$d1_tx90, na.rm = TRUE))
p2 <- rnorm(length(CW$d2_tx90), mean(CW$d2_tx90, na.rm = TRUE), sd(CW$d2_tx90, na.rm = TRUE))
p3 <- rnorm(length(CW$d1_tn90), mean(CW$d1_tn90, na.rm = TRUE), sd(CW$d1_tn90, na.rm = TRUE))
p4 <- rnorm(length(CW$d2_tn90), mean(CW$d2_tn90, na.rm = TRUE), sd(CW$d2_tn90, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_CW_sumtx90p.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Days with Tx >90p", ylim=c(0,15), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p2, prob=FALSE, breaks=5, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1,15, "Centralwestern region", cex=0.9, col="grey20", adj = 0)
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
dev.off()
graphics.off()
quartz()
png(filename = paste("rnorm_CW_sumtn90p.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Days with Tn >90p", ylim=c(0,15), main=NULL, xlim=c(0,1200), ylab = "Stations")
hist (p4, prob=FALSE, breaks=5, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("topright", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1,15, "Centralwestern region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()