Merge pull request #21 from felipe-treistman/validacoes

Validacoes
ONSBR · Feb 22, 2024 · 5d7d009 · 5d7d009
2 parents a5b7857 + d5e3de7
commit 5d7d009
Show file tree

Hide file tree

Showing 7 changed files with 48 additions and 12 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: smapOnsR
 Title: SMAP/ONS
-Version: 0.2.0
+Version: 0.2.1
 Authors@R: 
     person("Felipe", "Treistman", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-9948-8680"))

diff --git a/NEWS.md b/NEWS.md
@@ -1,8 +1,20 @@
+# smapOnsR 0.2.1
+
+## Major changes
+
+* Ajuste da precipitacao observada apos o processo de assimilacao limitado ao valor maximo de 1
+
+## Minor changes
+
+* escrita de arquivos de otimizacao e funcao_objetivo na funcao de execucao
+
+* ajuste nas funcoes de leitura de arquivo de saida
+
 # smapOnsR 0.2.0
 
 ## Major changes
 
-* Possibilidade de ajuste da precipitacao observada apos o processo de assimilacao. O processo de assimilacao permanece igual. Apos o termino da assimilacao a precipitação observada diaria é corrigida pelo fator do penultimo dia de assimilacao.
+* Possibilidade de ajuste da precipitacao observada apos o processo de assimilacao. O processo de assimilacao permanece igual. Apos o termino da assimilacao a precipitação observada diaria e corrigida pelo fator do penultimo dia de assimilacao.
 
 ## Minor changes
 

diff --git a/R/assimilacao.R b/R/assimilacao.R
@@ -135,7 +135,7 @@ assimilacao_oficial <- function(vetor_modelo, area, EbInic, TuInic, Supin, preci
     vetor_inicializacao <- unlist(inicializacao)
 
     if (ajusta_precipitacao == 1) {
-      ajuste$par[numero_dias_assimilacao] <- ajuste$par[(numero_dias_assimilacao - 1)]
+      ajuste$par[numero_dias_assimilacao] <- min(ajuste$par[(numero_dias_assimilacao - 1)], 1)
     }
 
     precipitacao_ponderada <- precipitacao_ponderada * ajuste$par[1:numero_dias_assimilacao]
@@ -150,7 +150,7 @@ assimilacao_oficial <- function(vetor_modelo, area, EbInic, TuInic, Supin, preci
     colnames(otimizacao) <- "otimizacao"
     otimizacao[, limite_inferior := limite_inferior]
     otimizacao[, limite_superior := limite_superior]
-    otimizacao[, variavel := c(paste0("prec (t-", 31:1,")"), "Ebin", "Supin")]
+    otimizacao[, variavel := c(paste0("prec (t-", 31:1, ")"), "Ebin", "Supin")]
 
     saida <- list(ajuste = ajuste, simulacao = simulacao, otimizacao = otimizacao)
     saida
@@ -384,7 +384,7 @@ assimilacao_evapotranspiracao <- function(vetor_modelo, area, EbInic, TuInic, Su
     vetor_inicializacao <- unlist(inicializacao)
 
     if (ajusta_precipitacao == 1) {
-      ajuste$par[numero_dias_assimilacao] <- ajuste$par[(numero_dias_assimilacao - 1)]
+      ajuste$par[numero_dias_assimilacao] <- min(ajuste$par[(numero_dias_assimilacao - 1)], 1)
     }
 
     precipitacao_ponderada <- precipitacao_ponderada * ajuste$par[1:numero_dias_assimilacao]

diff --git a/R/execucao.R b/R/execucao.R
@@ -24,11 +24,16 @@ executa_caso_oficial <- function(pasta_caso){
     entrada <- le_arq_entrada(pasta_entrada)
 
     saida <- rodada_encadeada_oficial(entrada$parametros,
-        entrada$inicializacao, entrada$precipitacao, entrada$previsao_precipitacao, entrada$evapotranspiracao, entrada$vazao,
+        entrada$inicializacao, entrada$precipitacao, entrada$previsao_precipitacao, 
+        entrada$evapotranspiracao, entrada$vazao,
         entrada$postos_plu, entrada$datas_rodadas, entrada$caso$nome_subbacia)
 
     escreve_previsao(pasta_saida, saida$previsao)
     escreve_ajuste(pasta_saida, saida$assimilacao)
+    write.table(saida$otimizacao, file = file.path(pasta_saida, "otimizacao.csv"),
+                row.names = FALSE, quote = FALSE, sep = ";")
+    write.table(saida$funcao_objetivo, file = file.path(pasta_saida, "funcao_objetivo.csv"),
+                row.names = FALSE, quote = FALSE, sep = ";")
 }
 
 #' Executa caso oficial com entrada nova

diff --git a/R/leitura_saida.R b/R/leitura_saida.R
@@ -134,15 +134,16 @@ le_previsao_2 <- function(pasta_saida, data) {
     dat <- fread(arq)
     dat <- melt(dat, id.vars = 1:2, variable.name = "data_previsao", value.name = "vazao")
     dat[, data_previsao := as.Date(data_previsao, format = "%d/%m/%Y")]
-    dat[, data_caso := data]
+    dat[, data_caso := as.Date(data, format = "%d/%m/%Y")]
     dat[, variavel := "Qcalc"]
     colnames(dat) <- tolower(colnames(dat))
     setcolorder(dat, c("data_caso", "data_previsao", "modelo", "subbacia" , "variavel", "vazao"))
     colnames(dat)[3:4] <- c("cenario", "nome")
     colnames(dat)[6] <- c("valor")
     dat[, nome := tolower(nome)]
     dat[, cenario := tolower(cenario)]
-    setorder(dat, nome, data_caso, cenario, data_previsao)
+    data.table::setorder(dat, nome, data_caso, cenario, data_previsao)
+
 
     return(dat)
 }
@@ -174,7 +175,8 @@ le_previsao <- function(pasta_saida, cenario, sub_bacia, data) {
     previsao[, cenario := cenario]
     previsao[, variavel := "Qcalc"]
     previsao[, nome := sub_bacia]
+    previsao[, Data := as.Date(Data, format = "%d/%m/%Y")]
     data.table::setnames(previsao, c("Data", "Qcal"), c("data_previsao", "valor"))    
-
+    data.table::setcolorder(previsao, c("data_caso", "data_previsao", "nome", "cenario", "variavel", "valor"))
     return(previsao)
 }
diff --git a/R/rodadas_encadeadas.R b/R/rodadas_encadeadas.R
@@ -220,7 +220,8 @@ rodada_encadeada_oficial <- function(parametros, inicializacao, precipitacao_obs
             }
 
             if (ajusta_precipitacao == 1) {
-                precipitacao[data_previsao <= dataRodada, valor := valor * ajuste$ajuste$par[numero_dias_assimilacao - 1]]
+                precipitacao[data_previsao <= dataRodada, valor := 
+                valor * min(ajuste$ajuste$par[numero_dias_assimilacao - 1], 1)]
             }
             matriz_precipitacao_previsao <- array(precipitacao[data_previsao < (dataRodada + numero_dias_previsao + kt_max) & data_rodada == dataRodada & 
                 data_previsao >= (dataRodada - kt_min - 1), valor], c(numero_dias_previsao + kt_max + kt_min + 1, numero_cenarios))
@@ -502,7 +503,8 @@ rodada_encadeada_etp <- function(parametros, inicializacao, precipitacao_observa
             }
 
             if (ajusta_precipitacao == 1) {
-                precipitacao[data_previsao <= dataRodada, valor := valor * ajuste$ajuste$par[numero_dias_assimilacao - 1]]
+                precipitacao[data_previsao <= dataRodada, valor := 
+                valor * min(ajuste$ajuste$par[numero_dias_assimilacao - 1], 1)]
             }
             matriz_precipitacao_previsao <- array(precipitacao[data_previsao < (dataRodada + numero_dias_previsao + kt_max) & data_rodada == dataRodada & 
                 data_previsao >= (dataRodada - kt_min - 1), valor], c(numero_dias_previsao + kt_max + kt_min + 1, numero_cenarios))

diff --git a/tests/testthat/test-rodadas_encadeadas.R b/tests/testthat/test-rodadas_encadeadas.R
@@ -74,6 +74,21 @@ test_that("testa rodada ecmwf", {
   } else {
     expect_true(abs(round(saida$previsao[nome == "pimentalt" & variavel == "Qcalc" & cenario == "ecmwf_1", valor][31], 0) - 10166) < 10166 * 0.01)
   }
+
+  entrada$inicializacao[variavel == "ajusta_precipitacao", valor := 1]
+
+  set.seed(129852)
+  saida <- rodada_encadeada_oficial(entrada$parametros,
+      entrada$inicializacao, entrada$precipitacao, entrada$previsao_precipitacao, entrada$evapotranspiracao, entrada$vazao,
+      entrada$postos_plu, entrada$datas_rodadas, entrada$caso$nome_subbacia)
+
+  secao <- sessionInfo()
+
+  if (secao$R.version$os == "mingw32") {
+    expect_equal(round(saida$previsao[nome == "tucurui" & variavel == "Qcalc" & cenario == "ecmwf_1", valor][1], 0), 2211)
+  } else {
+    expect_true(abs(round(saida$previsao[nome == "tucurui" & variavel == "Qcalc" & cenario == "ecmwf_1", valor][1], 0) - 2211) < 2211 * 0.01)
+  }
 })
 
 test_that("testa rodada ecmwf formato oficial", {
@@ -146,7 +161,7 @@ test_that("testa rodada com serie temporal etp", {
     datas_rodadas$numero_dias_previsao <- datas_rodadas$numero_dias_previsao + 2
 
     inicializacao <- data.table::data.table(nome = c(rep("avermelha", 9), rep("ssimao2", 9)),
-    variavel = rep(c("Ebin", "Supin", "Tuin", "numero_dias_assimilacao", 
+    variavel = rep(c("Ebin", "Supin", "Tuin", "numero_dias_assimilacao",
     "limite_inferior_ebin", "limite_superior_ebin", "limite_superior_prec", "limite_inferior_prec",
     "ajusta_precipitacao"), 2),
     valor = c(218.71, 46.69, 0.2891, 31, 0.8, 1.2, 2, 0.5, 0, 441.67, 256.98, 0.3141, 31, 0.8, 1.2, 2, 0.5, 0))