30-gamma.R

# =================================================================
# 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 (Exemplo para Janeiro) ---
dados <- read_delim("C:/Users/Jefferson B. Bezerra/Downloads/ETo/dez.csv", 
                    delim = ";", escape_double = FALSE, trim_ws = TRUE)

## Rows: 930 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 - OPÇÃO A (Mês Fixo 30 Dias)
# =================================================================
agrupar_eto_fixo <- function(df, tamanho_janela) {
  df %>%
    group_by(ano) %>%
    mutate(dia_mes = day(data)) %>%
    # DESCARTA O DIA 31 PARA PADRONIZAR
    filter(dia_mes <= 30) %>% 
    mutate(grupo = ceiling(dia_mes / tamanho_janela)) %>% 
    group_by(ano, grupo) %>%
    summarise(ETo_sum = sum(Eto, na.rm = TRUE), .groups = "drop")
}

# Novos intervalos coerentes com 30 dias
nomes_fixos <- c("30 Dias", "15 Dias", "10 Dias", "6 Dias", "5 Dias", "3 Dias")

df_lista <- list(
  "30 Dias" = agrupar_eto_fixo(dados, 30),
  "15 Dias" = agrupar_eto_fixo(dados, 15),
  "10 Dias" = agrupar_eto_fixo(dados, 10),
  "6 Dias"  = agrupar_eto_fixo(dados, 6),
  "5 Dias"  = agrupar_eto_fixo(dados, 5),
  "3 Dias"  = agrupar_eto_fixo(dados, 3)
)

# Unificar dados
df_bruto_total <- bind_rows(df_lista, .id = "Escala")
df_bruto_total$Escala <- factor(df_bruto_total$Escala, levels = nomes_fixos)

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

# 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 (HISTOGRAMA + GAMA + ESTATÍSTICAS)
# =================================================================

grafico_validacao <- 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.5, family = "mono", fontface = "bold") +
  facet_wrap(~Escala, scales = "free", ncol = 2) + 
  scale_fill_viridis_d(option = "turbo") +
  labs(
    title = "Ajuste da Distribuição Gama (Mês Padronizado: 30 dias)",
    subtitle = "Escalas Temporais com Intervalos Homogêneos",
    x = "Evapotranspiração de Referência (mm)",
    y = "Densidade"
  ) +
  theme_bw(base_size = 11) +
  theme(
    legend.position = "none",
    strip.background = element_rect(fill = "gray20", color = "black"),
    strip.text = element_text(face = "bold", size = 10, color = "white"),
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold", hjust = 0.5, size = 14)
  )

# =================================================================
# 5. EXPORTAÇÃO
# =================================================================
ggsave("ETo_Gama_30dias_dez.png", plot = grafico_validacao, type = "cairo", 
       width = 20, height = 24, units = "cm", dpi = 600)

## Warning: Using ragg device as default. Ignoring `type` and `antialias`
## arguments

print(grafico_validacao)

30-gamma.R

Jefferson B. Bezerra

2026-02-12