-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathANA2Swat.Rmd
147 lines (130 loc) · 7.59 KB
/
ANA2Swat.Rmd
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
---
title: "ANAdatabasetoSwat"
author: "Lima, M. Edberto"
date: "29/04/2021"
output: html_document
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
library(RODBC) ### Acesss database conection
library(tidyr)
library(dplyr)
library(lubridate) ### Packpage to convert data
library(tidyverse) ###
library(plyr) ###
library(plotly) ###
library(ggplot2) ###
library(zoo)
```
```{r Precipitation}
setwd("C:/temp/SwatProject_Grb") ### Select folder
# connect to Access Data base file
db <- ".\\Gauges\\ANA\\BcoFluPlu.mdb" ### Database files (.mdb) can be download form https://www.snirh.gov.br/hidroweb/ or requested by email - hidro@ana.gov.br
con <- odbcConnectAccess2007(db)
### Select weather gauge Stations
qryPCP <- "SELECT * FROM Chuvas WHERE (year(Data) >= 2008) AND
(EstacaoCodigo = 2851022 OR
EstacaoCodigo = 2951010 OR
EstacaoCodigo = 2951070 OR
EstacaoCodigo = 2851044 OR
EstacaoCodigo = 2851021 OR
EstacaoCodigo = 2851024 OR
EstacaoCodigo = 2852046 OR
EstacaoCodigo = 2852050 OR
EstacaoCodigo = 2952001 OR
EstacaoCodigo = 2952038 OR
EstacaoCodigo = 2851052 OR
EstacaoCodigo = 2852052 OR
EstacaoCodigo = 2852053 )
order by EstacaoCodigo , Data"
### Get access to precipation table
TblPrecipitation <- sqlQuery(con, qryPCP)
### Format and reorder the preciptation data
lstPrecipitation <- lapply(1:nrow(TblPrecipitation), function(x){
output_df <- data.frame("Code" = TblPrecipitation[x,"EstacaoCodigo"],
"NivelConsistencia" = TblPrecipitation[x,"NivelConsistencia"],
"Date" = seq.Date(floor_date(as.Date(TblPrecipitation[x,"Data"], "%d/%m/%Y", tz = '' ),"month"),
by = "day", length.out = days_in_month(as.Date(TblPrecipitation[x,"Data"], "%d/%m/%Y", tz = ''))),
"TipoMedicaoChuvas" = TblPrecipitation[x,"TipoMedicaoChuvas"],
"PCP" = round(t(TblPrecipitation[x,] %>%
select(starts_with("Chuva"),-contains("Status")) %>%
select(1:days_in_month(as.Date(TblPrecipitation[x,"Data"], "%d/%m/%Y", tz = '')))),2),
"ChuvaStatus" = (t(TblPrecipitation[x,] %>%
select(starts_with("Chuva") & ends_with("Status")) %>%
select(1:days_in_month(as.Date(TblPrecipitation[x,"Data"], "%d/%m/%Y", tz = ''))))))
names(output_df)[5:6] <- c("PCP", "ChuvaStatus")
return(output_df)
})
### Unlist
PCP.ANA <- do.call(rbind, lstPrecipitation)
### Order PCP by gauge stations and fill date gaps
List.WeatherStaGrb <- PCP.ANA %>%
group_split(Code)
### Fill gaps and format lists to produce a Heatmap graph
List.WeatherStaGrb <- lapply(1:length(List.WeatherStaGrb), function(x){
List.WeatherStaGrb[[x]] %>%
complete(Date = seq(ymd("2008-01-01"), #### Fill the gaps on data
ymd("2020-12-31"),1),
fill = list(Code = as.character(unique(List.WeatherStaGrb[[x]][,"Code"])),
PCP = -99)) %>% ### Riquired by SWAT
mutate( year = factor(year(Date)), ### To creat a calendar organize months and weeks
weekday = wday(Date, label = T, week_start = 1),
yearmonth = factor(as.yearmon(Date)),
month = month(Date, label = T, abbr = F),
wk = as.POSIXlt(Date)$wday,
week = as.numeric(isoweek(Date)),
day = day(Date)) %>%
mutate(week = case_when(month == "December" & week == 1 ~ 53, ### Correct weeks of year
month == "January" & week %in% 52:53 ~ 0,
TRUE ~ week),
pcat = cut(PCP, c(-1, 0, .5, 1:5, 7, 9, 15, 20, 25, 30, 200)))%>% ### Create class of PCP
ddply(.(yearmonth),transform,monthweek=1+week-min(week)) ### Convert weeks of year to month
})
### Set color to Heatmap
pubu <- RColorBrewer::brewer.pal(9, "BuPu")[3:9]
col_p <- colorRampPalette(pubu)
### Create a heatmapCalendar
lapply(1:length(List.WeatherStaGrb), function(x){
ggplot(List.WeatherStaGrb[[x]],
aes(weekday, monthweek, fill = pcat)) +
geom_tile(colour = "white", size = .4) +
guides(fill = guide_colorsteps(barwidth = 15,
barheight = .3,
title.position = "top")) +
scale_fill_manual(values = c("gray90", col_p(13)),
na.value = "grey65", drop = FALSE) +
facet_grid(year~ month, scales = "free_y") +
labs(subtitle = paste("Weather Station code:",unique(List.WeatherStaGrb[[x]]$Code)),
title = "Daily precipitation",
caption = "Source: ANA - Agência Nacional de Águas e Saneamento Básico",
fill = "mm") +
scale_y_reverse() +
theme_classic()+
theme(legend.position = "bottom",
legend.text = element_text( hjust = .5),
legend.title = element_text(size = 8, hjust = 1),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 90,hjust=1, vjust=0.2),
panel.grid = element_blank(),
panel.background = element_blank(),
plot.caption = element_text( hjust = 1, size = 8),
text = element_text(size=8))
# ggsave(file= paste0(unique(dat_pr[[x]]$Code),"_PCP.png"), ### Save files
# width = 26, height = 18, units = "cm",dpi = 300)
})
### Write a SwatFile - Need Check the initial date
lapply(1:length(List.WeatherStaGrb), function(x){
write.table(List.WeatherStaGrb[[x]][,c("PCP")],
col.names = 20080501, ### Start date
row.names = F,
quote = F,
sep = ",",
file=paste0("PCP_",unique(List.WeatherStaGrb[[x]][1,"Code"]),"_",
format(as.Date(List.WeatherStaGrb[[x]][1,"Date"]),"%Y%m%d"),"_",
format(as.Date(List.WeatherStaGrb[[x]][dim(List.WeatherStaGrb[[x]])[1],"Date"],"Date"),"%Y%m%d"),
".txt"))
})
```