-
Notifications
You must be signed in to change notification settings - Fork 0
/
hw_rnorm_hist_Tx_slopes_explore.R
228 lines (168 loc) · 10.4 KB
/
hw_rnorm_hist_Tx_slopes_explore.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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
setwd("//Volumes/Extended/INMET/allstations/")
dfm <- read.csv("slopes.csv", header = TRUE, sep = ",")
summary(dfm)
attach(dfm)
par(mfrow=c(2,2))
hist(d2_tx_slope)
boxplot(d2_tx_slope, d1_tx_slope)
graphics.off()
boxplot(dfm[,2:9])
plot(d1_tx_slope, d1_tn_slope)
plot(d1_txs_slope, d1_tns_slope)
qqnorm(dfm[,2])
qqline(dfm[,2])
boxplot(dfm[,2])
ba <- summary(d1_tn_slope)
250*(ba["Median"])
p1 <- rnorm(length(dfm$d1_tx_slope), mean(dfm$d1_tx_slope, na.rm = TRUE), sd(dfm$d1_tx_slope, na.rm = TRUE))
p2 <- rnorm(length(dfm$d2_tx_slope), mean(dfm$d2_tx_slope, na.rm = TRUE), sd(dfm$d2_tx_slope, na.rm = TRUE))
p3 <- rnorm(length(dfm$d1_tn_slope), mean(dfm$d1_tn_slope, na.rm = TRUE), sd(dfm$d1_tn_slope, na.rm = TRUE))
p4 <- rnorm(length(dfm$d2_tn_slope), mean(dfm$d2_tn_slope, na.rm = TRUE), sd(dfm$d2_tn_slope, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_All_txslope.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Slope of Tx vs. time", ylim=c(0,100), main=NULL, xlim=c(-1.5,1), 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")
col_code <- c(rgb(0,0,0.5,0.4), rgb(0.5,0,0,0.4))
legend("left", 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_tnslope.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Slope of Tn vs. time", ylim=c(0,100), main=NULL, xlim=c(-1.5,1), 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("left", 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_tx_slope), mean(NE$d1_tx_slope, na.rm = TRUE), sd(NE$d1_tx_slope, na.rm = TRUE))
p2 <- rnorm(length(NE$d2_tx_slope), mean(NE$d2_tx_slope, na.rm = TRUE), sd(NE$d2_tx_slope, na.rm = TRUE))
p3 <- rnorm(length(NE$d1_tn_slope), mean(NE$d1_tn_slope, na.rm = TRUE), sd(NE$d1_tn_slope, na.rm = TRUE))
p4 <- rnorm(length(NE$d2_tn_slope), mean(NE$d2_tn_slope, na.rm = TRUE), sd(NE$d2_tn_slope, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_NE_txslope.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Slope of Tx vs. time", ylim=c(0,50), main=NULL, xlim=c(-1.5,1), ylab = "Stations")
hist (p2, prob=FALSE, breaks=3, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1.5,50, "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("left", 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_tnslope.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Slope of Tn vs. time", ylim=c(0,50), main=NULL, xlim=c(-1.5,1), 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("left", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1.5,50, "Northeast region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- North -----------------
p1 <- rnorm(length(N$d1_tx_slope), mean(N$d1_tx_slope, na.rm = TRUE), sd(N$d1_tx_slope, na.rm = TRUE))
p2 <- rnorm(length(N$d2_tx_slope), mean(N$d2_tx_slope, na.rm = TRUE), sd(N$d2_tx_slope, na.rm = TRUE))
p3 <- rnorm(length(N$d1_tn_slope), mean(N$d1_tn_slope, na.rm = TRUE), sd(N$d1_tn_slope, na.rm = TRUE))
p4 <- rnorm(length(N$d2_tn_slope), mean(N$d2_tn_slope, na.rm = TRUE), sd(N$d2_tn_slope, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_N_txslope.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=3, col=rgb(0.5,0,0,0.4), xlab="Slope of Tx vs. time", ylim=c(0,20), main=NULL, xlim=c(-1.5,1), 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.5,20, "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("left", 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_tnslope.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Slope of Tn vs. time", ylim=c(0,20), main=NULL, xlim=c(-1.5,1), ylab = "Stations")
hist (p4, prob=FALSE, breaks=3, 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("left", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1.5,20, "North region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- Southeast -----------------
p1 <- rnorm(length(SE$d1_tx_slope), mean(SE$d1_tx_slope, na.rm = TRUE), sd(SE$d1_tx_slope, na.rm = TRUE))
p2 <- rnorm(length(SE$d2_tx_slope), mean(SE$d2_tx_slope, na.rm = TRUE), sd(SE$d2_tx_slope, na.rm = TRUE))
p3 <- rnorm(length(SE$d1_tn_slope), mean(SE$d1_tn_slope, na.rm = TRUE), sd(SE$d1_tn_slope, na.rm = TRUE))
p4 <- rnorm(length(SE$d2_tn_slope), mean(SE$d2_tn_slope, na.rm = TRUE), sd(SE$d2_tn_slope, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_SE_txslope.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Slope of Tx vs. time", ylim=c(0,40), main=NULL, xlim=c(-1.5,1), ylab = "Stations")
hist (p2, prob=FALSE, breaks=3, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1.5,40, "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("left", 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_tnslope.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Slope of Tn vs. time", ylim=c(0,40), main=NULL, xlim=c(-1.5,1), 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("left", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1.5,40, "Southeast region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- South -----------------
p1 <- rnorm(length(S$d1_tx_slope), mean(S$d1_tx_slope, na.rm = TRUE), sd(S$d1_tx_slope, na.rm = TRUE))
p2 <- rnorm(length(S$d2_tx_slope), mean(S$d2_tx_slope, na.rm = TRUE), sd(S$d2_tx_slope, na.rm = TRUE))
p3 <- rnorm(length(S$d1_tn_slope), mean(S$d1_tn_slope, na.rm = TRUE), sd(S$d1_tn_slope, na.rm = TRUE))
p4 <- rnorm(length(S$d2_tn_slope), mean(S$d2_tn_slope, na.rm = TRUE), sd(S$d2_tn_slope, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_S_txslope.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Slope of Tx vs. time", ylim=c(0,20), main=NULL, xlim=c(-1.5,1), 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.5,20, "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("left", 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_tnslope.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=5, col=rgb(0.5,0,0,0.4), xlab="Slope of Tn vs. time", ylim=c(0,20), main=NULL, xlim=c(-1.5,1), ylab = "Stations")
hist (p4, prob=FALSE, breaks=3, 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("left", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1.5,20, "South region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()
#--------------- Centralwestern -----------------
p1 <- rnorm(length(CW$d1_tx_slope), mean(CW$d1_tx_slope, na.rm = TRUE), sd(CW$d1_tx_slope, na.rm = TRUE))
p2 <- rnorm(length(CW$d2_tx_slope), mean(CW$d2_tx_slope, na.rm = TRUE), sd(CW$d2_tx_slope, na.rm = TRUE))
p3 <- rnorm(length(CW$d1_tn_slope), mean(CW$d1_tn_slope, na.rm = TRUE), sd(CW$d1_tn_slope, na.rm = TRUE))
p4 <- rnorm(length(CW$d2_tn_slope), mean(CW$d2_tn_slope, na.rm = TRUE), sd(CW$d2_tn_slope, na.rm = TRUE))
quartz()
png(filename = paste("rnorm_CW_txslope.png"), width = 500, height= 500, res = 90)
hist (p1, prob=FALSE, breaks=10, col=rgb(0.5,0,0,0.4), xlab="Slope of Tx vs. time", ylim=c(0,15), main=NULL, xlim=c(-1.5,1), ylab = "Stations")
hist (p2, prob=FALSE, breaks=2.5, col=rgb(0,0,0.5,0.4), add=TRUE)
leg_txt <- c("1997-2006","2007-2016")
text(-1.5,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("left", 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_tnslope.png"), width = 500, height= 500, res = 90)
hist (p3, prob=FALSE, breaks=6, col=rgb(0.5,0,0,0.4), xlab="Slope of Tn vs. time", ylim=c(0,15), main=NULL, xlim=c(-1.5,1), ylab = "Stations")
hist (p4, prob=FALSE, breaks=6, 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("left", legend = leg_txt, fill=col_code, horiz = FALSE, cex = 1, border="black", box.col = NULL, bty = "n")
text(-1.5,15, "Centralwestern region", cex=0.9, col="grey20", adj = 0)
dev.off()
graphics.off()