# Carregando funções
source("SCRIPTS/funcoes_roli.R")
source("~/Dropbox/Dissertacao/scripts/auxs_funs.R")
# Carregando dados
data_obs <-
read_rds("../Dados/Dados_SM_Nov2013Set2015_1hora_posFBE.rds") %>%
as.tbl %>%
select(date,Tar_fill,UR_fill,Rg_fill,Li,Le) %>%
dplyr::rename(Ta = Tar_fill,Rg = Rg_fill, rh = UR_fill) %>%
mutate(Rpot = Rg_Rpot(date)) %>%
mutate(es = (rh/100) * (0.6112 * exp((17.67*Ta)/(243.5+Ta))) *10) %>%
mutate(Ta = Ta+273.15)
## Loading required package: REddyProc
## Warning in Rg_Rpot(structure(c(1384905600, 1384909200, 1384912800,
## 1384916400, : Latitude e Longitude de Santa Maria
data_obs$K <-
kloudines(Rg = data_obs %>%
select(date,Rg) %>%
dplyr::rename(Rg = Rg)) %>%
select(K) %>%
unlist %>%
as.vector()
## Warning in kloudines(Rg = data_obs %>% select(date, Rg) %>%
## dplyr::rename(Rg = Rg)): Latitude e Longitude de Santa Maria
timePlot(data_obs %>% mutate(daytime = to.daylight(date,lon=-53.76,lat=-29.72,timezone=-4)),
"K", type = "daytime")
sigma <- 5.669 * 10e-9
cloud_par <- c("CQB","CKC","CCB","CKZ","CWU","CJG")
nuv_par <- c("FAN","FBR","FHY","FKZ","FIJ")
aten_par <- c("EAN","EBR","ESW","EIJ","EBT","EID","EKZ","EPR","ABM","ALH","AGB")
roli_comb <- rbind(expand.grid("-",aten_par), expand.grid(cloud_par,nuv_par) )
all_params <-
lapply(1:nrow(roli_comb), function(i){
# i = 3
pars_ <- roli_comb[i,] %>% unlist %>% as.vector
roli_i(data = data_obs,E_fun = pars_[2],C_fun = pars_[1])
}) %>% bind_cols()
tdy.roli <-
cbind(data_obs %>% select(date,Li,Ta,es),all_params) %>%
gather(params,value,c(-date,-Li, -Ta,-es)) %>%
as.tbl
tdy.roli%>%
select(date,Li,params,value) %>%
timePlot(c("Li","value"),type = "params", layout = c(4,11))
Seleção
Erro do instrumento
sel_params <-
tdy.roli %>%
group_by(params, dia = as.Date(date)) %>%
dplyr::summarise(obs = mean(Li, na.rm = TRUE), sim = mean(value, na.rm = TRUE)) %>%
ungroup() %>%
group_by(params) %>%
dplyr::summarise(obs_avg = mean(obs, na.rm = TRUE),
RMSEd = rmse( sim = sim, obs = obs, na.rm = TRUE),
pbias = pbias(sim = sim, obs = obs, na.rm = TRUE)
) %>%
filter(RMSEd <= obs_avg*0.1 & pbias < 5 & pbias > -5 ) %>%
select(params) %>%
t() %>% c()
sel_params
## [1] "AGB_-" "EID_-" "FAN_CQB" "FBR_CCB" "FBR_CJG" "FBR_CKZ" "FBR_CQB"
## [8] "FBR_CWU" "FHY_CCB" "FHY_CJG" "FHY_CKC" "FHY_CKZ" "FHY_CQB" "FHY_CWU"
## [15] "FIJ_CCB" "FIJ_CJG" "FIJ_CQB" "FIJ_CWU" "FKZ_CKC" "FKZ_CKZ"
tdy.dly.roli <-
tdy.roli %>%
select(date,Li, params,value) %>%
filter(params %in% sel_params) %T>%
timePlot(c("Li","value") , type = "params", lty = 1, col = c("red","gray50"), layout = c(4,5)) %>%
group_by(params, dia = as.Date(date)) %>%
dplyr::summarise(obs = mean(Li), sim = mean(value)) %>%
ungroup() %>% dplyr::rename(date=dia) %T>%
timePlot(c("obs","sim") , type = "params",lty = 1, col = c("red","gray25"), layout = c(4,5))
taylor_all <-
TaylorDiagram(tdy.dly.roli %>% as.data.frame(),
obs = "obs",mod = "sim",
group = "params",key.columns = 3,
pch = 5:10, col = 2:7, cex = 1,
main = "Série média diária da ROLi")
Atenuação
Céu claro
tdy.clr.roli <-
tdy.roli %>%
filter(params %in% sel_params) %>%
spread(params,value) %>%
mutate(Rg = data_obs$Rg %>% as.numeric, Rpot = data_obs$Rpot %>% as.numeric) %>%
gather(params,sim,c(-date,-Li,-Rg,-Rpot,-Ta,-es)) %>%
dplyr::rename(obs = Li) %>%
filter(Rpot > 0.1) %>%
filter(Rg > 0.75*Rpot) %T>%
timePlot(c("Rg","Rpot"),
plot.type = "h",
y.relation = "same",
yscale.components = yscale.components.subticks)
TaylorDiagram(tdy.clr.roli %>% as.data.frame(),
obs = "obs",mod = "sim",
group = "params",key.columns = 1,
pch = 5:10, col = 2:7, cex = 1,
main = "Série horários de céu claro ROLi")
Céu nublado
tdy.ocst.roli <-
tdy.roli %>%
filter(params %in% sel_params) %>%
spread(params,value) %>%
mutate(Rg = data_obs$Rg %>% as.numeric, Rpot = data_obs$Rpot %>% as.numeric) %>%
gather(params,sim,c(-date,-Li,-Rg,-Rpot,-Ta,-es)) %>%
dplyr::rename(obs = Li) %>%
filter(Rpot > 0.1) %>%
filter(Rg <= 0.75*Rpot) %T>%
timePlot(c("Rg","Rpot"),
plot.type = "h",
y.relation = "same",
yscale.components = yscale.components.subticks)
TaylorDiagram(tdy.ocst.roli %>% as.data.frame(),
obs = "obs",mod = "sim",
group = "params",key.columns = 1,
pch = 5:10, col = 2:7, cex = 1,
main = "Série horários de céu nublado ROLi")
Estação do ano
# add column with estção do ano
tdy.sly.roli <-
tdy.roli %>%
filter(params %in% sel_params) %>%
mutate(daytime = to.daylight(date,lon =-53.76, lat = -29.72,timezone = -4)) %>%
cutData(type = "season",hemisphere = "southern") %>%
spread(params,value) %>%
mutate(Rg = data_obs$Rg %>% as.numeric, Rpot = data_obs$Rpot %>% as.numeric) %>%
gather(params,sim,c(-date,-Li,-Rg,-Rpot,-Ta,-es,-season,-daytime)) %>%
dplyr::rename(obs = Li)
TaylorDiagram(tdy.sly.roli %>% as.data.frame(),
obs = "obs",mod = "sim",
group = "params",type = c("daytime","season"),
key.columns = 1,
pch = 5:10, col = 2:7, cex = 1,
main = "Série ...",layout = c(2,4))
##
## Using 'season' in data frame and not date-based openair version.
## This may result in different behaviour compared with openair calculations.
Decomposição do MSE
# Avaliação das parametrizações::
stats.roli <- eval.params(data_ = tdy.roli)
# NSE
# stats.roli.NSE <-
select_stats(roli_list = stats.roli,idx = "NSE")
# RMSE
# stats.roli.RMSE <-
select_stats(roli_list = stats.roli,idx = "RMSE")
# PBIAS
# stats.roli.PBIAS <-
select_stats(roli_list = stats.roli,idx = "PBIAS %")
## Estatística da série:
estats.roli <-
lapply(unique(tdy.roli$params), function(i){ # i = "FBM_CQB"
tdy.roli.filt <-
tdy.roli %>%
filter(params == i)
gof.data <-
gof(sim = tdy.roli.filt$value,obs = tdy.roli.filt$Li,na.rm = TRUE) %>%
as.data.frame() %>%
set_names(i)
}) %>% bind_cols()
estats.roli$stats <- rownames(gof(1:10,10:1))
estats.roli %<>%
gather(params,value,-stats) %>%
arrange(stats)
estats.roli %>%
separate(params,sep = "_",into = c("emis","aten")) %>%
spread(emis,value) %>%
filter(stats == "MAE")
tdy.roli %>%
select(date,params,value) %>%
group_by(day = day(date),params) %>%
mutate_each(funs(mean),c(-date,-params)) %>%
ungroup %>%
ggplot( aes(x=date,y=value,colour = params)) +
geom_line()
# NOTAS :
#
sel_params <-
spread.tdy %>%
group_by(params, dia = as.Date(date)) %>%
dplyr::summarise(obs = mean(Li, na.rm = TRUE), sim = mean(value, na.rm = TRUE)) %>%
ungroup() %>%
group_by(params) %>%
dplyr::summarise(obs_avg = mean(obs, na.rm = TRUE), RMSEd = rmse(sim = sim, obs = obs, na.rm = TRUE)) %>%
filter(RMSEd <= obs_avg*0.1) %>%
select(params) %>%
t() %>% c()
#
# FAZER TAYLOR TaylorDiagram() das series estimadas
spread.tdy <-
tdy.roli %>%
select(date,Li,params,value) %>%
# spread(params,value) %>%
separate(params, into = c("emis","cloud"), sep = "_") %>%
filter( (substr(emis,1,1) == "E" & substr(cloud,1,2) == "CQ") | substr(emis,1,1) == "F") %>%
unite(col = params,c(emis,cloud))
taylor_all <- TaylorDiagram(spread.tdy %>% as.data.frame(),
obs = "Li",mod = "value",
group = "params",key.columns = 3, pch = 1:25, cex = 1)
new.update <-
spread.tdy %>%
group_by(params) %>%
dplyr::summarise(pbias = hydroGOF::pbias(obs = Li,sim= value), RMSE = rmse(obs = Li,sim= value)) %>%
dplyr::rename(MyGroupVar = params) %>%
full_join(taylor_all$data, by = "MyGroupVar") %>%
mutate(MyGroupVar = factor(MyGroupVar))
taylor_all$data <- new.update
d <- data.frame(taylor_all$data)
scatterPlot(d, x = "pbias", y = "RMSE", cols = sample(colors(), 34) )
# criterio R2 > 0.5, explica mais de 50% da variabilidade da Li obs
by_rmse <- d %>%
filter(MyGroupVar %in% sel_params) %>%
filter(R**2 >= 0.5 & (pbias < 5 & pbias > -5)) %>%
#arrange(pbias) %>%
arrange(RMSE)
by_pbias <- d %>%
filter(MyGroupVar %in% sel_params) %>%
filter(R**2 >= 0.5 & (pbias < 5 & pbias > -5)) %>%
arrange(abs(pbias))
by_rmse[with(by_rmse, order(abs(pbias), RMSE)), ]
by_rmse %>%
mutate(pbias_abs = abs(pbias)) %>%
doBy::orderBy(~ RMSE + pbias_abs,data = .)