Estimativa da Radiação de Onda Longa incidente Paraíso do Sul.

# Carregando funções
source("~/Dropbox/NEW_WORKS/PPGMET_2016/Rmd/SCRIPTS/funcoes_roli.R")
source("~/Dropbox/Dissertacao/scripts/auxs_funs.R")

# Carregando dados
data_obs <- 
  read_rds("../Dados/radiation_data_HL_PRS.rds") %>%
    mutate(Rpot = Rg_Rpot(date,lon = -53.18,lat = -29.71,timezone = -4))  %>%
    mutate(Ta = Temp_C+ 273.15,
           rh = RH,
           Rg = SWin,
           Li = LWin,
           es = (rh/100) * (0.6112 * exp((17.67*Temp_C)/(243.5+Temp_C))) *10) %>%
    select(date,Ta,rh,es,Rpot,Rg, Li) %T>%
  timePlot(c("Rpot","Rg","Ta","rh","es","Li"))

## Loading required package: REddyProc

## Warning: attributes are not identical across measure variables; they will
## be dropped

# 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)

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

data_obs %>% 
    mutate(daytime = to.daylight(date,lon = -53.18,lat = -29.71,timezone=-4)) %>%
    timePlot("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] "EAN_-"   "EBT_-"   "EID_-"   "EIJ_-"   "EPR_-"   "FKZ_CCB" "FKZ_CJG"
##  [8] "FKZ_CKC" "FKZ_CKZ" "FKZ_CQB" "FKZ_CWU"

tdy.dly.roli <- 
    tdy.roli %>%
    select(date,Li, params,value) %T>%
    # filter(params %in% sel_params) %T>% 
    timePlot(c("Li","value") , type = "params", lty = 1, col = c("red","gray50"), layout = c(3,4)) %>%
    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(3,4)) 

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

days.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) %>%
    group_by(day = as.Date(date)) %>%
    dplyr::summarize(R = cor(x = sim,y=obs)) %>%
    filter(R > 0.4) %>%
    ungroup() %>%
    select(day) %>%
    t %>% c

tdy.dly.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( format(date,"%Y-%m-%d") %in% days.clr.roli)  %>%
    filter(Rpot > 0.1 ) %T>%
    timePlot(c("Rg","Rpot"),
             plot.type = "h",
             y.relation = "same", 
             yscale.components = yscale.components.subticks) 

TaylorDiagram(tdy.dly.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

days.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) %>%
    group_by(day = as.Date(date)) %>%
    dplyr::summarize(R = cor(x = sim,y=obs)) %>%
    filter(R < 0.4) %>%
    ungroup() %>%
    select(day) %>%
    t %>% c

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( format(date,"%Y-%m-%d") %in% days.ocst.roli)  %T>%
    timePlot(c("Rpot","Rg"),plot.type = "h",
             col = c("black","red"), lty = 1,
             y.relation = "same", group = TRUE,
             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

MSE_nd <- function(data_in){
    sim = data_in$sim
    obs = data_in$obs
    idx <- valindex(sim,obs)
    nidx <- length(idx)
    
    Dmean  <-  (mean(sim,na.rm = TRUE) - mean(obs,na.rm = TRUE))^2
    Dsd <- (sd(sim,na.rm = TRUE) - sd(obs,na.rm = TRUE))^2
    return(c(Dmean = Dmean,Dsd = Dsd))
}

data_in <- 
tdy.roli %>%
    group_by(params) %>%
    select(Li,value) %>%
    rename(obs = Li,sim = value) 

 # 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 = .)