# =================================================================
# 1. CARREGAMENTO E CONFIGURAÇÕES
# =================================================================
library(tidyverse)
## Warning: pacote 'tidyverse' foi compilado no R versão 4.4.3
## Warning: pacote 'ggplot2' foi compilado no R versão 4.4.3
## Warning: pacote 'tidyr' foi compilado no R versão 4.4.1
## Warning: pacote 'purrr' foi compilado no R versão 4.4.1
## Warning: pacote 'dplyr' foi compilado no R versão 4.4.1
## Warning: pacote 'stringr' foi compilado no R versão 4.4.1
## Warning: pacote 'forcats' foi compilado no R versão 4.4.1
## Warning: pacote 'lubridate' foi compilado no R versão 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(viridis)
## Warning: pacote 'viridis' foi compilado no R versão 4.4.3
## Carregando pacotes exigidos: viridisLite
## Warning: pacote 'viridisLite' foi compilado no R versão 4.4.1
library(Cairo)
## Warning: pacote 'Cairo' foi compilado no R versão 4.4.3
# --- Importação (Arquivo de Fevereiro) ---
dados <- read_delim("C:/Users/Jefferson B. Bezerra/Downloads/ETo/fev.csv",
delim = ";", escape_double = FALSE, trim_ws = TRUE)
## Rows: 847 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (1): data
## dbl (6): Tmax, Tmin, UR, n, U2, Eto
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dados <- dados %>% mutate(data = dmy(data), ano = year(data))
# =================================================================
# 2. FUNÇÃO DE AGRUPAMENTO - FEVEREIRO (Fixo 28 Dias)
# =================================================================
agrupar_eto_feb <- function(df, tamanho_janela) {
df %>%
group_by(ano) %>%
mutate(dia_mes = day(data)) %>%
# CONSIDERA APENAS ATÉ O DIA 28 (IGNORA DIA 29 EM ANOS BISSEXTOS)
filter(dia_mes <= 28) %>%
mutate(grupo = ceiling(dia_mes / tamanho_janela)) %>%
group_by(ano, grupo) %>%
summarise(ETo_sum = sum(Eto, na.rm = TRUE), .groups = "drop")
}
# Escalas perfeitamente divisíveis por 28
nomes_feb = c("28 Dias", "14 Dias", "7 Dias", "4 Dias")
df_lista <- list(
"28 Dias" = agrupar_eto_feb(dados, 28),
"14 Dias" = agrupar_eto_feb(dados, 14),
"7 Dias" = agrupar_eto_feb(dados, 7),
"4 Dias" = agrupar_eto_feb(dados, 4)
)
# Unificar dados
df_bruto_total <- bind_rows(df_lista, .id = "Escala")
df_bruto_total$Escala <- factor(df_bruto_total$Escala, levels = nomes_feb)
# =================================================================
# 3. CÁLCULOS E ESTATÍSTICAS
# =================================================================
calc_estatisticas <- function(df_escala, nome) {
vetor <- df_escala$ETo_sum
m_ar <- mean(vetor, na.rm = TRUE)
A <- log(m_ar) - mean(log(vetor), na.rm = TRUE)
alfa <- (1 + sqrt(1 + (4/3) * A)) / (4 * A)
beta <- m_ar / alfa
ks <- ks.test(vetor, "pgamma", shape = alfa, scale = beta)
txt <- sprintf("α: %.2f\nβ: %.2f\nD: %.3f\np: %.3f",
alfa, beta, ks$statistic, ks$p.value)
return(data.frame(Escala = nome, Alfa = alfa, Beta = beta, Rotulo = txt))
}
resumo_estatistico <- map2_df(df_lista, names(df_lista), calc_estatisticas)
## Warning in ks.test.default(vetor, "pgamma", shape = alfa, scale = beta): não
## devem existir empates no teste de Kolmogorov-Smirnov de apenas uma amostra
## Warning in ks.test.default(vetor, "pgamma", shape = alfa, scale = beta): não
## devem existir empates no teste de Kolmogorov-Smirnov de apenas uma amostra
## Warning in ks.test.default(vetor, "pgamma", shape = alfa, scale = beta): não
## devem existir empates no teste de Kolmogorov-Smirnov de apenas uma amostra
resumo_estatistico$Escala <- factor(resumo_estatistico$Escala, levels = nomes_feb)
# Curvas teóricas
df_curvas <- resumo_estatistico %>%
group_by(Escala) %>%
do({
vetor_ref <- df_lista[[.$Escala]]$ETo_sum
x_seq <- seq(0, max(vetor_ref) * 1.1, length.out = 300)
data.frame(x = x_seq, densidade = dgamma(x_seq, shape = .$Alfa, scale = .$Beta))
}) %>%
ungroup()
# =================================================================
# 4. GRÁFICO FINAL (FEVEREIRO)
# =================================================================
grafico_feb <- ggplot() +
geom_histogram(data = df_bruto_total, aes(x = ETo_sum, y = after_stat(density), fill = Escala),
bins = 15, color = "white", alpha = 0.5) +
geom_line(data = df_curvas, aes(x = x, y = densidade), color = "black", linewidth = 0.8) +
geom_text(data = resumo_estatistico, aes(label = Rotulo),
x = Inf, y = Inf, hjust = 1.1, vjust = 1.2,
size = 3.8, family = "mono", fontface = "bold") +
facet_wrap(~Escala, scales = "free", ncol = 2) +
scale_fill_viridis_d(option = "turbo") +
labs(
title = "Ajuste da Distribuição Gama - Fevereiro (Fixo 28 dias)",
subtitle = "Intervalos Homogêneos: 28, 14, 7 e 4 dias",
x = "Evapotranspiração de Referência (mm)",
y = "Densidade"
) +
theme_bw(base_size = 12) +
theme(
legend.position = "none",
strip.background = element_rect(fill = "gray20", color = "black"),
strip.text = element_text(face = "bold", size = 11, color = "white"),
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5, size = 14)
)
# =================================================================
# 5. EXPORTAÇÃO
# =================================================================
ggsave("ETo_Gama_28dias_fev.png", plot = grafico_feb, type = "cairo",
width = 20, height = 20, units = "cm", dpi = 600)
## Warning: Using ragg device as default. Ignoring `type` and `antialias`
## arguments
print(grafico_feb)
