diff --git a/DESCRIPTION b/DESCRIPTION index 19d17fb..bb56743 100644 --- 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", , "felipe.treistman@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9948-8680")) diff --git a/NEWS.md b/NEWS.md index 9653d2c..a538394 100644 --- 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 index ade3688..1cbd6dd 100644 --- 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 index dc8bcf3..fdcf2e2 100644 --- 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 index 153edde..91bf466 100644 --- a/R/leitura_saida.R +++ b/R/leitura_saida.R @@ -134,7 +134,7 @@ 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")) @@ -142,7 +142,8 @@ le_previsao_2 <- function(pasta_saida, data) { 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 index ccc5cc1..4044818 100644 --- 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 index 8bb4dd6..3d78526 100644 --- 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))