-
Notifications
You must be signed in to change notification settings - Fork 1
/
PCP_Q
272 lines (223 loc) · 13 KB
/
PCP_Q
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
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
# Instruction to Download Monitoring Data for Canadian watersheds:
#1.register at GitHub under university login, indicate a student status to receive free membership.
#2. R libraries to download data from Hydat and Weather Canada are stored at GitHub.
install.packages("devtools") # If not already installed
devtools::install_github("ropensci/weathercan")
devtools::install_github("ropensci/tidyhydat")
#or
devtools::install_github("paleolimbot/hydatr")
#3. Using ArcGIS select nearby stations for flow and weather;
HydatStn<-c("02HM002", "02HM007", "O2HM004","02HM001","02HM010", "02HL102","02HM003", "02HE003", "02HE001", "02HM011", "02HM006", "02HM005", "02KF020", "02LA02302LA017", "02KF016")
#4. Use R libraries above to download data for those stations,
#Important - data should be used as data.frames, so use conversion statement :
x<-data.frame(x)
#5. Prepare a column for dates 1970-1-1 till today - it will server as a primary ID for connect with downloaded data from different time periods .
#code
TimeAll=as.integer(365.25*49) #we create sequence of time data for 49 years in daily scale from 1970s till Dec. 31, 2018
df_timeAll <- data.frame(date1 = seq(as.Date('1970-01-01'),by = 'day', len =TimeAll))
#all missing values will receive NA, which we shall change later for NAN for Winbugs use.
#5. Using primary ID column you can merge it with downloaded data from each downloaded station.
# Also rename columne to include station ID if necessary.
# merge two data frames by date column "ID"
total <- merge(df_timeAll,df_Hydat02EC005,by="PK",all.x = TRUE)
#6. Save the dataframe as an object for use in Winbugs :
dput(total, file = "PCP_Q", control = c("keepNA", "keepInteger", "niceNames", "showAttributes"))
#
dget(file, keep.source = FALSE)
# downloading flows
HydatStn <- c("02HM002", "02HM007","02HM001","02HM010", "02HL102","02HM003", "02HE003", "02HE001", "02HM011", "02HM006", "02HM005", "02LA017", "02KF016")
total <- data.frame(date1 = seq(as.Date('1970-01-01'),by = 'day', len =TimeAll))
names(total)[1] <- 'Date'
add_flow <- function(stnNo,total)
{
temp_q <- tidyhydat::hy_daily_flows(station_number = stnNo, start_date = "1970-01-01", end_date = "2018-12-31")
names(temp_q)[4] <- stnNo
temp_q$STATION_NUMBER <- NULL
temp_q$Parameter <- NULL
temp_q$Symbol <- NULL
total <- merge(total,temp_q,by = "Date",all.x = TRUE)
return(total)
}
for (i in 1:length(HydatStn)) {
print(HydatStn[i])
stnNo <- HydatStn[i]
total <- add_flow(stnNo,total) }
dput(total,file="total.dat")
#Assess the amount of missing values:
naData<-matrix(0, nrow=3, ncol=ncol(total))
for (i in 1:NCOL(total)) {
naData[1,i]<-nrow(total)
naData[2,i]<-sum(is.na(total[,i]))
naData[3,i]<-(naData[1,i]-naData[2,i])/(naData[1,i]}
#rounding the number of existing values as ration (1=100%).
round(naData[3,],2)
#######
#Extration of weather:
install.packages("weathercan")
library("weathercan")
#weatherstn<-c("TRENTON A", "BELLEVILLE", "BLOOMFIELD", "PICTON", "CRESSY", "KINGSTON A", "CATARAQUI TS", "HARTINGTON", "GODFREY")
weatherstn<-c(5126, 4859, 4860, 4861, 4862, 4863, 5005, 5006, 4911, 7671, 4295, 52985, 50428, 4242, 4287, 4284) #unique EC station ids
TimeAll=as.integer(365.25*49) #we create sequence of time data for 49 years in daily scale from 1970s till Dec. 31, 2018
df_timeAll <- data.frame(date1 = seq(as.Date('1970-01-01'),by = 'day', len =TimeAll))
Final <-data.frame(date1 = seq(as.Date('1970-01-01'),by = 'day', len =TimeAll))
add_weather <- function(stnNa,final)
+ {
+ kam_q <- weather_dl(station_ids = stnNa, start_date = "1970-01-01", end_date = "2018-12-31")
+ names(kam_q)[4] <- stnNa
+ kam_q$station_name <-NULL
+ kam_q&Parameter <-NULL
+ kam_q&Symbol <-NULL
+ final <-merge (final, kam_q, by = "Date", all.x=TRUE)
+ return (final)
+ }
#weather data # extract a list of weather stations
StationID<- rclimateca::ec_climate_locations_all
#plot stations by stationID
write.csv(list,file = "weatherStn.csv")
#Extract all weather stations:
StationID <- rclimateca::ec_climate_locations_all
list <- data.frame(BoQ_Weather_Stn)
stn <- list$station_id
kam_q <- weathercan::weather_dl(station_ids = stn[1], interval = "day", start = "1970-01-01", end = "2018-12-31")
kam_q <- data.frame(kam_q)
datePcpc <- data.frame(kam_q[,11])
datePcpc$PCP <- kam_q$total_precip
===== Dec.26 2018 ====
getwd()
total<-dget(file="total.dat", keep.source = FALSE)
#extract stationID (internal ID of Weather Canada)
StationID <- rclimateca::ec_climate_locations_all
write.csv(StationID,file = "StationID_Clima.csv")
#plot in ArcGIS for further analysis and clipping stations of interest
#read a table with stations of interest into R
BoQ_Weather_Stn <- read.table(file = "BoQ_Weather_Stn.txt")
list <- data.frame(BoQ_Weather_Stn)
weatherstn<-c(5126,30723, 4859, 4860, 4861, 4862,4863, 5005, 5006, 4911,7671, 4295, 50428,52985, 4242, 4287, 4284) #unique EC station ids
Stn_names<- c("TRENTON A", "TRENTON AWOS", "BELLEVILLE", "BELLEVILLE OWRC", "BELLEVILLE PAR LAB","BLOOMFIELD", "BLOOMFIELD WEST", "PICTON", "PICTON A", "CRESSY", "KINGSTON (AUT)", "KINGSTON A1", "KINGSTON A2", "KINGSTON A3", "CATARAQUI TS", "HARTINGTON", "GODFREY")
names(weatherstn) <- c(Stn_names)
#Adding pcp as a function
add_pcp <- function(stnNo,total)
{
temp_q <- weather_dl(station_ids = stnNo, start = "1970-01-01", end = "2018-12-24", interval="day", string_as = NULL)
#temp_q <-temp_q[,c("station_name", "station_id", "date", "total_precip", "total_rain", "total_snow")]
temp_q <-temp_q[,c("date", "total_precip")]
names(temp_q) <- c("Date",stnNo)
#names(temp_q) <- c("Date",paste0("tot_pcp_",stnNo))
total <- merge(total,temp_q,by = "Date",all.x = TRUE)
return(total)
}
#actual processing a list of stations
for (i in 1:length(weatherstn)) {
print(weatherstn[i])
stnNo <- weatherstn[i]
total <- add_pcp(stnNo,total) }
################## Dec 27, 2018
dates <- total$Date
#.Define a new variable(dates) by using date column
water_year <- function(dates, start_month=10) {
# Define a function (water_year) from October to September
dates.posix = as.POSIXlt(dates)
# Convert dates into POSIXlt
offset = ifelse(dates.posix$mon >= start_month - 1, 1, 0)
# Year offset
adj.year = dates.posix$year + 1900 + offset
# Water year
adj.year
}
# Return the water year
wateryear <- water_year(dates,10)
#specify that oct marks the first month of the hydrological year, and sept therefore marks the last month.
total = cbind(total, wateryear)
#add the hydrological year as a column in the df
Month<-format(total$Date,"%m")
#Extract month from the date column in the df
total= cbind(total, Month)
#Add Month to the df
getwd()
QPCP <- data.frame(dget(file = "total2.dat"))
library(lubridate)
QPCP$Month <- lubridate::month(QPCP$Date)
QPCP$Year <- lubridate::year(QPCP$Date)
#QPCP$YearWater <- 0
library(dplyr)
# QPCP$YearWater <- QPCP %>% mutate(Year_water = ifelse(.$Month %in% 6:12, .$Year, .$Year - 1))
#
#rm(QPCP)
# adding water year /hydrological year
QPCP <- QPCP %>% mutate(Year_water = ifelse(.$Month %in% 1:6, .$Year + 1, .$Year))
#adding growth period : 1- yes, 0 - no.
QPCP <- QPCP %>% mutate(PeriodGrowth = ifelse(.$Month %in% 6:10, 1, 0))
QPCP %>%
group_by(.$Year_water)
summarise_all(funs(sum(!is.na(.)))) -> qpcp.totCount #-> qpcp.totCount
QPCP %>% group_by(.$Year_water) %>% summarise_all(funs(sum(!is.na(.)))) -> pivot.CountWaterYear
QPCP %>% filter(.$PeriodGrowth == 1) %>% group_by(.$Year_water) %>% summarise_all(funs(sum(!is.na(.)))) -> pivot.CountGrowth
pivot.CountGrowth <- data.frame(pivot.CountGrowth)
pivot.Growth.Binary <- data.frame(pivot.CountGrowth)
pivot.Growth.Binary[pivot.Growth.Binary != 153] <- 0
pivot.Growth.Binary[pivot.Growth.Binary == 153] <- 1
pivot.Growth.Binary$Year_water <- pivot.CountGrowth$..Year_water
pivot.Growth.Binary <- data.frame(subset(pivot.Growth.Binary,pivot.Growth.Binary$PeriodGrowth == 1))
colNames <- colnames(pivot.Growth.Binary)
write.csv(colNames,file = "colNames.csv")
#pivot.Growth.Binary %>% select(.$Year_water,.$X02HM002.x,.$X5126) %>% filter(.$X02HM002.x == 1 & .$X5126 == 1)
library(dplyr)
pivot.Growth.Binary %>% select_(.dots = c("X02HM002.x","X5126","Year_water")) %>% filter(.$X02HM002.x == 1 & .$X5126 == 1) -> dateTemp
yearTemp <- dateTemp$Year_water
# now we need to filter data from those years in yearTemp
QPCP %>% filter(.$Year_water %in% dateTemp$Year_water) -> subsetTemp
subsetTemp$ln_q_X02HM002 <- log10(subsetTemp$X02HM002.x + 0.001)
library(zoo)
library(DataCombine)
subsetTemp = subsetTemp %>%
mutate(pcp_5126_2d = rollmean(x = .$X5126, 2, align = "right", fill = NA))
# mutate(pcp_5126_2d = rollmean(x = .$X5126, 2, align = "center", fill = NA))
subsetTemp$ln_pcp_5126_2d <- log10(subsetTemp$pcp_5126_2d + 0.001)
subsetTemp %>% select(ln_pcp_5126_2d, ln_q_X02HM002) %>% filter(.$ln_pcp_5126_2d !=-3) -> dataPlot
library(reshape2)
input.subsetTemp <- melt(dataPlot)
library(gcookbook)
sp <- ggplot(dataPlot, aes(x = dataPlot$ln_pcp_5126_2d , y = dataPlot$ln_q_X02HM002)) + geom_point(shape = 1, alpha = 1/4, position = position_jitter(width = 1, height = 1)) + geom_smooth(color = "red")
#stat_smooth(method = lm, se = FALSE, colour = "black")
sp
ggsave(sp, filename = "5126_2d vs 02HM002.png")
#to plot two figures on one layout use this package:
library("gridExtra")#
library(ggplot2)
library(gridExtra)
library(grid)
sp0<-ggplot(total3b, aes(x= total3b$`5126`, y = total3b$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp1<- ggplot(MA, aes(x= MA$`5126 2d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp2<- ggplot(MA, aes(x= MA$`5126 3d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp3<- ggplot(MA, aes(x= MA$`5126 4d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp4<- ggplot(MA, aes(x= MA$`5126 5d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp5<- ggplot(MA, aes(x= MA$`5126 6d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp6<- ggplot(MA, aes(x= MA$`5126 7d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp7<- ggplot(MA, aes(x= MA$`5126 8d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp8<- ggplot(MA, aes(x= MA$`5126 9d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp9<- ggplot(MA, aes(x= MA$`5126 10d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp10<- ggplot(MA, aes(x= MA$`5126 11d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp11<- ggplot(MA, aes(x= MA$`5126 12d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp12<- ggplot(MA, aes(x= MA$`5126 13d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp13<- ggplot(MA, aes(x= MA$`5126 14d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp14<- ggplot(MA, aes(x= MA$`5126 15d` , y = MA$`02HM002`)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp14a<- ggplot(total3b, aes(x= total3b$Ln5126, y = total3b$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp15 <-ggplot(MA, aes(x= MA$`Ln5126 2d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp16 <-ggplot(MA, aes(x= MA$`Ln5126 3d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp17 <-ggplot(MA, aes(x= MA$`Ln5126 4d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp18 <-ggplot(MA, aes(x= MA$`Ln5126 5d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp19 <-ggplot(MA, aes(x= MA$`Ln5126 6d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp20 <-ggplot(MA, aes(x= MA$`Ln5126 7d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp21 <-ggplot(MA, aes(x= MA$`Ln5126 8d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp22 <-ggplot(MA, aes(x= MA$`Ln5126 9d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp23 <-ggplot(MA, aes(x= MA$`Ln5126 10d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp24 <-ggplot(MA, aes(x= MA$`Ln5126 11d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp25 <-ggplot(MA, aes(x= MA$`Ln5126 12d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp26 <-ggplot(MA, aes(x= MA$`Ln5126 13d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp27 <-ggplot(MA, aes(x= MA$`Ln5126 14d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
sp28 <-ggplot(MA, aes(x= MA$`Ln5126 15d` , y = MA$Ln02HM002)) + geom_point(size=0.5) + geom_smooth(color = "red")
#construct scatterplots. Find and replace all to change the station names and codes.
panel <- grid.arrange(sp0,sp1,sp2,sp3,sp4,sp5,sp6,sp7,sp8,sp9,sp10,sp11,sp12,sp13,sp14,sp14a,sp15,sp16,sp17,sp18,sp19,sp20,sp21,sp22,sp23,sp24,sp25,sp26,sp27,sp28, ncol = 6, top="DEPOT CREEK AT BELLROCK (02HM002) vs. TRENTON A (5126)")
ggsave(panel, file = "5126 vs 02HM002 .png", width = 10, height = 8)
dev.off()
# Save plot