df_raw <- read_parquet(PATH_OCORR_CLIMA)

df_raw <- df_raw %>% left_join(segments, by="segmento")

df <- df_raw |>
  filter(coord_valida) |>
  mutate(
    ano     = year(data_inicio),
    mes     = month(data_inicio),
    mes_lab = month(data_inicio, label = TRUE, abbr = TRUE),
    evento_precip_ext  = tp_above_p95      == 1,
    evento_vento_forte = ws10_above_60kmh  == 1,
    evento_calor       = dia_calor_35      == 1,
    evento_gelo        = ciclo_gelo_degelo == 1,
    qualquer_extremo   = evento_precip_ext | evento_vento_forte |
                         evento_calor      | evento_gelo
  )

tipos_analisar <- df |>
  count(tipo_harmonizado, sort = TRUE) |>
  filter(n >= MIN_OCORR_TIPO) |>
  pull(tipo_harmonizado)

df_tipos <- df |> filter(tipo_harmonizado %in% tipos_analisar)

features_continuas <- c(
  "tp_d", "tp_h_max", "tp_acc_3h", "tp_acc_6h", "tp_acc_12h",
  "tp_acc_24h", "tp_acc_48h", "tp_acc_72h",
  "t2m_min", "t2m_max", "t2m_mean", "amplitude_termica",
  "ws10_max", "ws10_mean",
  "rh_min", "rh_max", "rh_mean",
  "indice_calor", "n_dias_seco_prev"
)

features_binarias <- c(
  "tp_above_p95", "tp_above_p99", "tp_above_20mm", "tp_h_above_10",
  "dia_calor_35", "dia_calor_40", "dia_frio", "ciclo_gelo_degelo",
  "onda_calor", "onda_frio", "ws10_above_60kmh", "ws10_above_80kmh",
  "ws10_above_100kmh", "rh_alta", "chuva_vento_forte", "precip_pos_seca"
)

1 Visão Geral do Dataset

liga si

tibble(
  Metrica = c(
    "Total de ocorrências",
    "Com coordenadas válidas",
    "Com features climáticas",
    "Período coberto",
    "Linhas distintas",
    "Segmentos distintos",
    "Tipos de ocorrência distintos",
    paste0("Tipos com >= ", MIN_OCORR_TIPO, " ocorrências (analisáveis)")
  ),
  Valor = c(
    format(nrow(df_raw), big.mark = " "),
    format(sum(df_raw$coord_valida, na.rm = TRUE), big.mark = " "),
    format(sum(!is.na(df_raw$tp_d)), big.mark = " "),
    paste(min(df$ano, na.rm = TRUE), "-", max(df$ano, na.rm = TRUE)),
    format(n_distinct(df$linha, na.rm = TRUE), big.mark = " "),
    format(n_distinct(df$segmento, na.rm = TRUE), big.mark = " "),
    format(n_distinct(df$tipo_harmonizado, na.rm = TRUE), big.mark = " "),
    as.character(length(tipos_analisar))
  )
) |>
  kable(align = c("l", "r"), booktabs = TRUE,
        col.names = c("Métrica", "Valor")) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    latex_options     = c("striped", "hold_position"),
    full_width        = FALSE
  )
Métrica Valor
Total de ocorrências 206 616
Com coordenadas válidas 204 593
Com features climáticas 204 593
Período coberto 2010 - 2025
Linhas distintas 61
Segmentos distintos 71
Tipos de ocorrência distintos 64
Tipos com >= 200 ocorrências (analisáveis) 37

2 Distribuição Base das Ocorrências

2.1 Evolução anual

df |>
  count(ano, grupo_harmonizado) |>
  ggplot(aes(ano, n, fill = grupo_harmonizado)) +
  geom_col() +
  scale_fill_manual(
    values = colorRampPalette(cores_principais)(n_distinct(df$grupo_harmonizado))
  ) +
  scale_x_continuous(breaks = seq(min(df$ano, na.rm = TRUE),
                                   max(df$ano, na.rm = TRUE), 2)) +
  scale_y_continuous(labels = comma) +
  labs(title = "Ocorrências por ano e grupo",
       x = NULL, y = "N.º de ocorrências", fill = "Grupo") +
  guides(fill = guide_legend(nrow = 3))

2.2 Volume por tipo de ocorrência

df |>
  count(tipo_harmonizado, grupo_harmonizado) |>
  mutate(tipo_harmonizado = fct_reorder(tipo_harmonizado, n)) |>
  ggplot(aes(n, tipo_harmonizado, fill = grupo_harmonizado)) +
  geom_col() +
  geom_vline(xintercept = MIN_OCORR_TIPO, linetype = "dashed",
             colour = "red", linewidth = 0.7) +
  annotate("text", x = MIN_OCORR_TIPO + 20, y = 3,
           label = paste0("Limiar análise (", MIN_OCORR_TIPO, ")"),
           colour = "red", size = 3, hjust = 0) +
  scale_fill_manual(
    values = colorRampPalette(cores_principais)(n_distinct(df$grupo_harmonizado))
  ) +
  scale_x_continuous(labels = comma) +
  labs(title    = "Volume de ocorrências por tipo",
       subtitle = "Linha vermelha: limiar mínimo para análise individual",
       x = "N.º de ocorrências", y = NULL, fill = "Grupo") +
  guides(fill = guide_legend(nrow = 3))

2.3 Sazonalidade mensal

df_tipos |>
  count(tipo_harmonizado, mes) |>
  group_by(tipo_harmonizado) |>
  mutate(pct = n / sum(n)) |>
  ungroup() |>
  mutate(tipo_harmonizado = fct_reorder(
    tipo_harmonizado,
    ifelse(mes %in% c(11, 12, 1, 2, 3), pct, 0),
    .fun = sum
  )) |>
  ggplot(aes(mes, tipo_harmonizado, fill = pct)) +
  geom_tile(colour = "white", linewidth = 0.3) +
  scale_x_continuous(
    breaks = 1:12,
    labels = c("Jan","Fev","Mar","Abr","Mai","Jun",
               "Jul","Ago","Set","Out","Nov","Dez")
  ) +
  scale_fill_gradient(low = "#EBF5FB", high = "#1F4E79",
                      labels = percent_format(accuracy = 1)) +
  labs(title    = "Distribuição mensal por tipo de ocorrência",
       subtitle = "Proporção do total anual",
       x = NULL, y = NULL, fill = "% do total") +
  theme(axis.text.x     = element_text(angle = 30, hjust = 1),
        axis.text.y     = element_text(size = 8),
        legend.position = "right")

2.4 Distribuição por segmento (top 20)

#  df |>
#   count(segmento, linha) |>
#   slice_max(n, n = 20) |>
#   mutate(segmento = fct_reorder(as.character(segmento), n)) %>%
#   mutate(n)

top_seg <- df |>
  group_by(segmento, linha) |>
  summarise(
    n = n(),
    extensao = first(extensao),
    .groups = "drop"
  ) |>
  mutate(n_por_ext = n / extensao) |>
  slice_max(n_por_ext, n = 20) |>
  mutate(segmento = fct_reorder(as.character(segmento), n_por_ext))

top_seg |>
  ggplot(aes(n_por_ext, segmento, fill = as.character(linha))) +
  geom_col() +
  scale_x_continuous(labels = comma) +
  scale_fill_manual(
    values = colorRampPalette(cores_principais)(n_distinct(top_seg$linha))
  ) +
  labs(title = "Top 20 segmentos com mais ocorrências relativas (ocorrências / extensão)",
       x = "N.º de ocorrências", y = "Segmento", fill = "Linha") +
  theme(legend.position = "right",
        axis.text.y     = element_text(size = 8))

3 Distribuição Geográfica das Ocorrências

# Fronteiras de Portugal (requer pacote rnaturalearth)
# install.packages(c("rnaturalearth", "rnaturalearthdata", "sf"))
portugal <- ne_countries(scale = "medium", country = "Portugal",
                          returnclass = "sf")

# Contagem de ocorrências por célula ERA5
ocorr_por_celula <- df |>
  filter(!is.na(era5_lat), !is.na(era5_lon)) |>
  count(era5_lat, era5_lon, name = "n_ocorrencias")

# Mapa base: ocorrências totais por célula ERA5
p_mapa_total <- ggplot() +
  geom_sf(data = portugal, fill = "#F5F5F5", colour = "#AAAAAA",
          linewidth = 0.5) +
  geom_point(
    data = ocorr_por_celula,
    aes(x = era5_lon, y = era5_lat,
        size  = n_ocorrencias,
        colour = n_ocorrencias),
    alpha = 0.85
  ) +
  scale_colour_gradient(low  = "#D6E4F0", high = "#1F4E79",
                        name = "N.º ocorrências",
                        labels = comma) +
  scale_size_continuous(range = c(1, 12),
                        name  = "N.º ocorrências",
                        labels = comma) +
  guides(colour = guide_legend(), size = guide_legend()) +
  coord_sf(xlim = c(-9.8, -6.0), ylim = c(36.8, 42.2)) +
  labs(
    title    = "Concentração de ocorrências por célula ERA5-Land",
    subtitle = "Cada ponto = 1 célula da grelha (~9 km) | tamanho e cor = n.º de ocorrências",
    x = "Longitude", y = "Latitude"
  ) +
  theme(
    panel.background = element_rect(fill = "#EBF5FB", colour = NA),
    legend.position  = "right"
  )

print(p_mapa_total)

# Mapa por grupo de ocorrência (facet)
ocorr_por_celula_grupo <- df |>
  filter(!is.na(era5_lat), !is.na(era5_lon)) |>
  count(era5_lat, era5_lon, grupo_harmonizado, name = "n_ocorrencias")

ggplot() +
  geom_sf(data = portugal, fill = "#F5F5F5", colour = "#AAAAAA",
          linewidth = 0.4) +
  geom_point(
    data = ocorr_por_celula_grupo,
    aes(x = era5_lon, y = era5_lat,
        size   = n_ocorrencias,
        colour = n_ocorrencias),
    alpha = 0.8
  ) +
  scale_colour_gradient(low  = "#D6E4F0", high = "#C00000",
                        name = "N.º ocorrências",
                        labels = comma) +
  scale_size_continuous(range = c(0.5, 8),
                        name  = "N.º ocorrências",
                        labels = comma) +
  guides(colour = guide_legend(), size = guide_legend()) +
  coord_sf(xlim = c(-9.8, -6.0), ylim = c(36.8, 42.2)) +
  facet_wrap(~ grupo_harmonizado, ncol = 3) +
  labs(
    title    = "Distribuição geográfica por grupo de ocorrência",
    subtitle = "Células ERA5-Land com pelo menos uma ocorrência do grupo",
    x = "Longitude", y = "Latitude"
  ) +
  theme(
    panel.background = element_rect(fill = "#EBF5FB", colour = NA),
    strip.text       = element_text(size = 7, face = "bold"),
    legend.position  = "bottom",
    axis.text        = element_text(size = 6)
  )

4 Sinal Climático: Dias com vs Sem Ocorrências

celulas_com_ocorr <- df |> distinct(era5_lat, era5_lon)

dias_com <- df |>
  filter(!is.na(tp_d)) |>
  distinct(era5_lat, era5_lon, data_dia) |>
  mutate(tem_ocorrencia = TRUE)

# clima_ref <- read_parquet(
#   PATH_CLIMA,
#   col_select = c("lat", "lon", "time",
#                  all_of(features_continuas),
#                  all_of(features_binarias))
# ) |>
#   rename(era5_lat = lat, era5_lon = lon, data_dia = time) |>
#   mutate(data_dia = as.Date(data_dia)) |>
#   semi_join(celulas_com_ocorr, by = c("era5_lat", "era5_lon")) |>
#   left_join(dias_com, by = c("era5_lat", "era5_lon", "data_dia")) |>
#   mutate(tem_ocorrencia = replace_na(tem_ocorrencia, FALSE))

clima_ref <- read_parquet(
  PATH_CLIMA,
  col_select = c("lat", "lon", "time",
                 all_of(features_continuas),
                 all_of(features_binarias))
) |>
  rename(era5_lat = lat, era5_lon = lon) |>
  mutate(data_dia = as.Date(time)) |>   # ← converter aqui, antes do join
  select(-time) |>
  semi_join(celulas_com_ocorr, by = c("era5_lat", "era5_lon")) |>
  left_join(dias_com, by = c("era5_lat", "era5_lon", "data_dia")) |>
  mutate(tem_ocorrencia = replace_na(tem_ocorrencia, FALSE))

cat(sprintf("Dias x celulas na referencia: %s\n",
            format(nrow(clima_ref), big.mark = " ")))
## Dias x celulas na referencia: 1 742 312
cat(sprintf("  Com ocorrencia : %s (%.1f%%)\n",
            format(sum( clima_ref$tem_ocorrencia), big.mark = " "),
            100 * mean( clima_ref$tem_ocorrencia)))
##   Com ocorrencia : 154 385 (8.9%)
cat(sprintf("  Sem ocorrencia : %s (%.1f%%)\n",
            format(sum(!clima_ref$tem_ocorrencia), big.mark = " "),
            100 * mean(!clima_ref$tem_ocorrencia)))
##   Sem ocorrencia : 1 587 927 (91.1%)

4.1 Distribuição de features climáticas: com vs sem ocorrência

features_plot <- c(
  "tp_d", "tp_acc_24h", "tp_acc_72h",
  "t2m_max", "ws10_max", "amplitude_termica",
  "n_dias_seco_prev", "rh_mean"
)

labs_features <- c(
  tp_d              = "Precipitação diária (mm)",
  tp_acc_24h        = "Precipitação acumulada 24h (mm)",
  tp_acc_72h        = "Precipitação acumulada 72h (mm)",
  t2m_max           = "Temperatura máxima (°C)",
  ws10_max          = "Velocidade máx. vento (m/s)",
  amplitude_termica = "Amplitude térmica (°C)",
  n_dias_seco_prev  = "Dias secos anteriores",
  rh_mean           = "Humidade relativa média (%)"
)

for (feat in features_plot) {
  cat(sprintf("\n\n### %s\n\n", labs_features[feat]))

  p <- clima_ref |>
    filter(
      !is.na(.data[[feat]]),
      .data[[feat]] < quantile(.data[[feat]], 0.99, na.rm = TRUE)
    ) |>
    mutate(grupo = if_else(tem_ocorrencia, "Com ocorrência", "Sem ocorrência")) |>
    ggplot(aes(x = grupo, y = .data[[feat]], fill = grupo)) +
    geom_violin(alpha = 0.7, draw_quantiles = c(0.25, 0.5, 0.75)) +
    geom_boxplot(width = 0.1, fill = "white", outlier.shape = NA) +
    scale_fill_manual(values = c(
      "Com ocorrência" = "#1F4E79",
      "Sem ocorrência" = "#AAAAAA"
    )) +
    labs(title = labs_features[feat], x = NULL, y = labs_features[feat]) +
    theme(legend.position = "none")

  print(p)
  cat("\n\n")
}

4.1.1 Precipitação diária (mm)

4.1.2 Precipitação acumulada 24h (mm)

4.1.3 Precipitação acumulada 72h (mm)

4.1.4 Temperatura máxima (°C)

4.1.5 Velocidade máx. vento (m/s)

4.1.6 Amplitude térmica (°C)

4.1.7 Dias secos anteriores

4.1.8 Humidade relativa média (%)

4.2 Elevação de risco (lift) por evento climático extremo

freq_base <- clima_ref |>
  summarise(across(all_of(features_binarias),
                   ~ mean(.x == 1, na.rm = TRUE))) |>
  pivot_longer(everything(), names_to = "feature", values_to = "freq_base")

freq_ocorr <- clima_ref |>
  filter(tem_ocorrencia) |>
  summarise(across(all_of(features_binarias),
                   ~ mean(.x == 1, na.rm = TRUE))) |>
  pivot_longer(everything(), names_to = "feature", values_to = "freq_ocorr")

risco <- freq_base |>
  left_join(freq_ocorr, by = "feature") |>
  mutate(
    lift        = freq_ocorr / freq_base,
    feature_lab = labs_binarias[feature]
  ) |>
  filter(!is.na(lift), !is.na(feature_lab))

risco |>
  mutate(feature_lab = fct_reorder(feature_lab, lift)) |>
  ggplot(aes(lift, feature_lab)) +
  geom_col(aes(fill = lift > 1.5), show.legend = FALSE) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40") +
  geom_text(aes(label = sprintf("%.2fx", lift)), hjust = -0.1, size = 3.2) +
  scale_fill_manual(values = c("TRUE" = "#1F4E79", "FALSE" = "#AAAAAA")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title    = "Elevação de risco (lift) por evento climático extremo",
    subtitle = "Rácio entre frequência do evento nos dias com ocorrência vs frequência base.\nValor > 1 indica sobre-representação de ocorrências nesses dias.",
    x = "Lift (rácio de elevação)", y = NULL
  )

5 Análise por Tipo de Ocorrência

5.1 Perfil climático por tipo

for (tipo in tipos_analisar) {

  cat(sprintf("\n\n### %s\n\n",
              str_trunc(tipo, width = 35, side = "right")))

  dias_tipo <- df |>
    filter(tipo_harmonizado == tipo, !is.na(tp_d)) |>
    distinct(era5_lat, era5_lon, data_dia) |>
    mutate(tem_tipo = TRUE)

  ref_tipo <- clima_ref |>
    left_join(dias_tipo, by = c("era5_lat", "era5_lon", "data_dia")) |>
    mutate(tem_tipo = replace_na(tem_tipo, FALSE))

  n_tipo <- sum(ref_tipo$tem_tipo)

  lift_tipo <- features_binarias |>
    map_dfr(function(feat) {
      fb <- mean(ref_tipo[[feat]] == 1, na.rm = TRUE)
      fo <- mean(ref_tipo[[feat]][ref_tipo$tem_tipo] == 1, na.rm = TRUE)
      tibble(
        feature     = feat,
        feature_lab = labs_binarias[feat],
        lift        = if_else(fb > 0, fo / fb, NA_real_)
      )
    }) |>
    filter(!is.na(lift), !is.na(feature_lab))

  titulo_plot <- str_wrap(sprintf("%s  (n=%d)", tipo, n_tipo), width = 45)

  p_lift <- lift_tipo |>
    mutate(feature_lab = fct_reorder(feature_lab, lift)) |>
    ggplot(aes(lift, feature_lab, fill = lift > 1.5)) +
    geom_col(show.legend = FALSE) +
    geom_vline(xintercept = 1, linetype = "dashed", colour = "grey40") +
    geom_text(aes(label = sprintf("%.2fx", lift)), hjust = -0.1, size = 2.5) +
    scale_fill_manual(values = c("TRUE" = "#1F4E79", "FALSE" = "#AAAAAA")) +
    scale_x_continuous(expand = expansion(mult = c(0, 0.2))) +
    labs(title = titulo_plot, subtitle = "Lift por evento extremo",
         x = "Lift", y = NULL) +
    theme(
      plot.title  = element_text(size = 8, face = "bold"),
      axis.text.y = element_text(size = 8)
    )

  p_tp <- ref_tipo |>
    filter(!is.na(tp_d),
           tp_d < quantile(tp_d, 0.99, na.rm = TRUE)) |>
    mutate(grupo = if_else(tem_tipo, "Com ocorrência", "Sem ocorrência")) |>
    ggplot(aes(tp_d, fill = grupo, colour = grupo)) +
    geom_density(alpha = 0.4) +
    scale_fill_manual(values = c(
      "Com ocorrência" = "#1F4E79",
      "Sem ocorrência" = "#AAAAAA"
    )) +
    scale_colour_manual(values = c(
      "Com ocorrência" = "#1F4E79",
      "Sem ocorrência" = "#555555"
    )) +
    labs(title  = "Precipitação diária (mm)",
         x = "tp_d (mm/dia)", y = "Densidade",
         fill = NULL, colour = NULL) +
    theme(legend.position  = "top",
          plot.title       = element_text(size = 9))

  print(p_lift + p_tp + plot_layout(widths = c(2, 1)))
  cat("\n\n")
}

5.1.1 INSTALAÇÕES DE SINALIZAÇÃO | CON…

5.1.2 INSTALAÇÕES DE SINALIZAÇÃO EM PA…

5.1.3 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

5.1.4 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

5.1.5 INSTALAÇÕES DE SINALIZAÇÃO | OCU…

5.1.6 RISCOS / ACIDENTES / INCIDENTES …

5.1.7 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

5.1.8 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

5.1.9 VIA | DEFICIÊNCIA DE VIA

5.1.10 IRREGULARIDADES NA EXECUÇÃO DE O…

5.1.11 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.12 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

5.1.13 PEDIDO DA EMPRESA FERROVIÁRIA

5.1.14 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.15 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.16 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.17 ESTRUTURAS | OCORRÊNCIA EM INSTA…

5.1.18 ERROS NOS PROCEDIMENTOS OPERACIO…

5.1.19 ORGANIZAÇÃO DO HORÁRIO TÉCNICO

5.1.20 OUTRAS CAUSAS

5.1.21 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.22 RISCOS / ACIDENTES / INCIDENTES …

5.1.23 VIA | AVARIA MECÂNICA EM AMV

5.1.24 INSTALAÇÕES DE ALIMENTAÇÃO ELÉCT…

5.1.25 INSTALAÇÕES DE ALIMENTAÇÃO ELÉCT…

5.1.26 INSTALAÇÕES DE TELECOMUNICAÇÕES

5.1.27 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.28 FALHA NOS SISTEMAS DE APOIO À EX…

5.1.29 INSTALAÇÕES DE ALIMENTAÇÃO ELÉCT…

5.1.30 INSTALAÇÕES DE SINALIZAÇÃO EM PA…

5.1.31 PESSOAL

5.1.32 INSTALAÇÕES DE ALIMENTAÇÃO ELÉCT…

5.1.33 OUTRAS CAUSAS | AVARIA EM MATERI…

5.1.34 ESTRUTURAS | DEFICIÊNCIA EM OUTR…

5.1.35 CONDIÇÕES METEREOLÓGICAS E CAUSA…

5.1.36 OBRAS PLANEADAS

5.1.37 FALHA NOS SISTEMAS DE APOIO À EX…

5.2 Impacto em dias climáticos extremos por tipo

df_tipos |>
  filter(!is.na(minutos_atraso), minutos_atraso > 0) |>
  mutate(
    condicao = case_when(
      evento_precip_ext  ~ "Precip extrema",
      evento_vento_forte ~ "Vento forte",
      evento_calor       ~ "Calor extremo",
      evento_gelo        ~ "Gelo-degelo",
      TRUE               ~ "Sem evento extremo"
    )
  ) |>
  ggplot(aes(x = condicao, y = minutos_atraso,
             fill = condicao, colour = condicao)) +
  geom_violin(alpha = 0.5, draw_quantiles = 0.5, show.legend = FALSE) +
  geom_boxplot(width = 0.15, fill = "white", outlier.shape = NA,
               show.legend = FALSE) +
  scale_y_log10(labels = comma) +
  scale_fill_manual(values = c(
    "Precip extrema"     = "#1F4E79",
    "Vento forte"        = "#2E75B6",
    "Calor extremo"      = "#ED7D31",
    "Gelo-degelo"        = "#70AD47",
    "Sem evento extremo" = "#AAAAAA"
  )) +
  scale_colour_manual(values = c(
    "Precip extrema"     = "#1F4E79",
    "Vento forte"        = "#2E75B6",
    "Calor extremo"      = "#ED7D31",
    "Gelo-degelo"        = "#70AD47",
    "Sem evento extremo" = "#555555"
  )) +
  facet_wrap(~ tipo_harmonizado, scales = "free_y", ncol = 4) +
  labs(
    title    = "Distribuição de minutos de atraso por condição climática e tipo",
    subtitle = "Escala logarítmica — linha central = mediana",
    x = NULL, y = "Minutos de atraso (log)"
  ) +
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1, size = 7),
    strip.text  = element_text(size = 7)
  )

6 Análise Espacial por Segmento

6.1 Concentração de ocorrências por segmento e tipo

top15_seg  <- df |> count(segmento) |>
  slice_max(n, n = 15) |> pull(segmento)
top15_tipo <- df |> count(tipo_harmonizado) |>
  slice_max(n, n = 15) |> pull(tipo_harmonizado)

df |>
  filter(segmento %in% top15_seg, tipo_harmonizado %in% top15_tipo) |>
  count(segmento, tipo_harmonizado) |>
  mutate(segmento = fct_reorder(as.character(segmento), n, .fun = sum)) |>
  ggplot(aes(tipo_harmonizado, segmento, fill = n)) +
  geom_tile(colour = "white", linewidth = 0.4) +
  geom_text(aes(label = n), size = 2.5, colour = "white") +
  scale_fill_gradient(low = "#D6E4F0", high = "#1F4E79") +
  scale_x_discrete(guide = guide_axis(angle = 40)) +
  labs(
    title = "Ocorrências por segmento x tipo (top 15 de cada)",
    x = NULL, y = "Segmento", fill = "N.º ocorrências"
  ) +
  theme(legend.position = "right",
        axis.text.x     = element_text(size = 7),
        axis.text.y     = element_text(size = 8))

6.2 Precipitação nas ocorrências por segmento

df |>
  filter(segmento %in% top15_seg, !is.na(tp_d)) |>
  mutate(segmento = as.character(segmento)) |>
  ggplot(aes(
    x    = tp_d,
    y    = fct_reorder(segmento, tp_d, .fun = median),
    fill = after_stat(x)
  )) +
  geom_density_ridges_gradient(scale = 2, rel_min_height = 0.01,
                                quantile_lines = TRUE, quantiles = 2) +
  scale_fill_gradient(low = "#EBF5FB", high = "#1F4E79") +
  scale_x_continuous(limits = c(0, 60)) +
  labs(
    title    = "Distribuição de precipitação diária nas ocorrências por segmento",
    subtitle = "Linha central = mediana | segmentos ordenados por mediana crescente",
    x = "Precipitação diária (mm)", y = "Segmento", fill = "mm"
  ) +
  theme(legend.position = "right",
        axis.text.y     = element_text(size = 8))

6.3 Curvas de sensibilidade ao limiar de precipitação

limiares <- c(1, 5, 10, 20, 30, 50)

for (tipo in tipos_analisar[1:min(5, length(tipos_analisar))]) {

  cat(sprintf("\n\n### %s\n\n",
              str_trunc(tipo, width = 35, side = "right")))

  segs_tipo <- df |>
    filter(tipo_harmonizado == tipo, !is.na(tp_d)) |>
    count(segmento) |>
    filter(n >= 20) |>
    pull(segmento)

  if (length(segs_tipo) < 3) {
    cat("*Segmentos insuficientes para análise (min. 3 com >= 20 ocorrências).*\n\n")
    next
  }

  # Calcular % acima de cada limiar — sem map_dfc para evitar erros de scope
  sens_limiares <- df |>
    filter(tipo_harmonizado == tipo,
           segmento %in% segs_tipo,
           !is.na(tp_d)) |>
    group_by(segmento) |>
    summarise(
      n_total        = n(),
      pct_acima_1mm  = mean(tp_d >= 1,  na.rm = TRUE),
      pct_acima_5mm  = mean(tp_d >= 5,  na.rm = TRUE),
      pct_acima_10mm = mean(tp_d >= 10, na.rm = TRUE),
      pct_acima_20mm = mean(tp_d >= 20, na.rm = TRUE),
      pct_acima_30mm = mean(tp_d >= 30, na.rm = TRUE),
      pct_acima_50mm = mean(tp_d >= 50, na.rm = TRUE),
      .groups = "drop"
    ) |>
    pivot_longer(
      starts_with("pct_acima"),
      names_to  = "limiar",
      values_to = "pct_ocorr_acima"
    ) |>
    mutate(
      limiar_mm = as.numeric(str_extract(limiar, "\\d+")),
      segmento  = as.character(segmento)
    )

  p <- sens_limiares |>
    ggplot(aes(limiar_mm, pct_ocorr_acima,
               colour = segmento, group = segmento)) +
    geom_line(linewidth = 1) +
    geom_point(size = 2.5) +
    scale_y_continuous(labels = percent_format(accuracy = 1),
                       limits = c(0, NA)) +
    scale_x_continuous(breaks = limiares) +
    scale_colour_manual(
      values = colorRampPalette(cores_principais)(length(segs_tipo))
    ) +
    labs(
      title    = str_wrap(sprintf("Curvas de sensibilidade — %s", tipo), 60),
      subtitle = "% de ocorrências do tipo neste segmento em dias com tp >= limiar.\nCurvas que caem mais lentamente = segmentos mais sensíveis.",
      x = "Limiar de precipitação (mm/dia)",
      y = "% ocorrências acima do limiar",
      colour = "Segmento"
    ) +
    theme(legend.position = "right",
          plot.title      = element_text(size = 10))

  print(p)
  cat("\n\n")
}

6.3.1 INSTALAÇÕES DE SINALIZAÇÃO | CON…

6.3.2 INSTALAÇÕES DE SINALIZAÇÃO EM PA…

6.3.3 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

6.3.4 INSTALAÇÕES DE SINALIZAÇÃO | AVA…

6.3.5 INSTALAÇÕES DE SINALIZAÇÃO | OCU…

7 Correlação entre Features Climáticas

mat_cor <- df |>
  select(all_of(features_continuas)) |>
  filter(complete.cases(across(everything()))) |>
  cor(method = "spearman")

rownames(mat_cor) <- colnames(mat_cor) <- recode(
  rownames(mat_cor),
  tp_d              = "Precip diária",
  tp_h_max          = "Precip max/h",
  tp_acc_3h         = "Acum 3h",
  tp_acc_6h         = "Acum 6h",
  tp_acc_12h        = "Acum 12h",
  tp_acc_24h        = "Acum 24h",
  tp_acc_48h        = "Acum 48h",
  tp_acc_72h        = "Acum 72h",
  t2m_min           = "T min",
  t2m_max           = "T max",
  t2m_mean          = "T media",
  amplitude_termica = "Amplitude T",
  ws10_max          = "Vento max",
  ws10_mean         = "Vento medio",
  rh_min            = "HR min",
  rh_max            = "HR max",
  rh_mean           = "HR media",
  indice_calor      = "Indice calor",
  n_dias_seco_prev  = "Dias secos ant."
)

corrplot(
  mat_cor,
  method      = "color",
  type        = "upper",
  order       = "hclust",
  col         = colorRampPalette(c("#1F4E79", "white", "#C00000"))(200),
  tl.cex      = 0.8,
  tl.col      = "black",
  addCoef.col = "black",
  number.cex  = 0.55,
  diag        = FALSE,
  title       = "Correlação de Spearman entre features climáticas",
  mar         = c(0, 0, 1.5, 0)
)

Nota de modelação: grupos de features com correlação > 0,8 são redundantes em modelos lineares. Para tree-based models a colinearidade é menos crítica mas dificulta a interpretação de importâncias.

8 Síntese e Implicações para Modelação

resumo_modelacao <- tipos_analisar |>
  map_dfr(function(tipo) {
    dias_tipo <- df |>
      filter(tipo_harmonizado == tipo, !is.na(tp_d)) |>
      distinct(era5_lat, era5_lon, data_dia) |>
      mutate(tem_tipo = TRUE)

    ref_t <- clima_ref |>
      left_join(dias_tipo, by = c("era5_lat", "era5_lon", "data_dia")) |>
      mutate(tem_tipo = replace_na(tem_tipo, FALSE))

    n_ocorr <- sum(ref_t$tem_tipo)
    n_segs  <- df |> filter(tipo_harmonizado == tipo) |>
      pull(segmento) |> n_distinct()

    top_feat <- features_binarias |>
      map_dfr(function(feat) {
        fb <- mean(ref_t[[feat]] == 1, na.rm = TRUE)
        fo <- mean(ref_t[[feat]][ref_t$tem_tipo] == 1, na.rm = TRUE)
        tibble(feature = feat,
               lift    = if_else(fb > 0, fo / fb, NA_real_))
      }) |>
      filter(!is.na(lift)) |>
      slice_max(lift, n = 1)

    tibble(
      tipo          = tipo,
      n_ocorrencias = n_ocorr,
      n_segmentos   = n_segs,
      feature_top   = labs_binarias[top_feat$feature],
      lift_top      = round(top_feat$lift, 2),
      modelavel     = n_ocorr >= MIN_OCORR_TIPO & top_feat$lift > 1.3
    )
  })

resumo_modelacao |>
  arrange(desc(lift_top)) |>
  mutate(modelavel = if_else(modelavel, "Sim", "")) |>
  kable(
    col.names = c("Tipo de Ocorrência", "N.º Ocorrências", "N.º Segmentos",
                  "Feature com maior lift", "Lift máx.", "Modelável"),
    align     = c("l", "r", "r", "l", "r", "c"),
    booktabs  = TRUE
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    latex_options     = c("striped", "hold_position", "scale_down"),
    full_width        = TRUE
  ) |>
  column_spec(6, bold = TRUE, color = "#1F4E79") |>
  row_spec(which(resumo_modelacao$lift_top > 1.5), background = "#EBF5FB")
Tipo de Ocorrência N.º Ocorrências N.º Segmentos Feature com maior lift Lift máx. Modelável
ESTRUTURAS &#124; OCORRÊNCIA EM INSTALAÇÕES 2817 48 Vento >= 60 km/h 35.68 Sim
ESTRUTURAS &#124; OCORRÊNCIA EM INSTALAÇÕES 2817 48 Chuva + vento forte 35.68 Sim
INSTALAÇÕES DE SINALIZAÇÃO EM PASSAGEM DE NÍVEL 483 19 Vento >= 60 km/h 34.69 Sim
INSTALAÇÕES DE SINALIZAÇÃO EM PASSAGEM DE NÍVEL 483 19 Chuva + vento forte 34.69 Sim
CONDIÇÕES METEREOLÓGICAS E CAUSAS NATURAIS &#124; TEMPESTADES (CHUVAS TORRENCIAIS/VENTOS) 232 18 Precip > P99 26.45 Sim
VIA &#124; AVARIA MECÂNICA EM AMV 1375 44 Vento >= 60 km/h 12.18 Sim
VIA &#124; AVARIA MECÂNICA EM AMV 1375 44 Chuva + vento forte 12.18 Sim
INSTALAÇÕES DE ALIMENTAÇÃO ELÉCTRICA &#124; FALTA DE ALIMENTAÇÃO PELA EDP 484 30 Precip > P99 5.46 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA EM MESA DE COMANDO 6932 45 Vento >= 60 km/h 4.83 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA EM MESA DE COMANDO 6932 45 Chuva + vento forte 4.83 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA ELÉTRICA EM AMV 13502 52 Temp max >= 40°C 3.73 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; OCUPAÇÃO INTEMPESTIVA EM CIRCUITO VIA 10937 53 Vento >= 60 km/h 3.06 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; OCUPAÇÃO INTEMPESTIVA EM CIRCUITO VIA 10937 53 Chuva + vento forte 3.06 Sim
VIA &#124; DEFICIÊNCIA DE VIA 5596 47 Temp max >= 40°C 3.04 Sim
INSTALAÇÕES DE ALIMENTAÇÃO ELÉCTRICA &#124; DEFICIÊNCIA DE CATENÁRIA 1007 45 Onda de calor 2.85 Sim
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; SDO(SISTEMA DE DETEÇÃO DE OBSTÁCULOS) 1799 27 Precip > P99 2.82 Sim
INSTALAÇÕES DE ALIMENTAÇÃO ELÉCTRICA &#124; DESARMES 684 21 Precip intensa pos-seca 2.68 Sim
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; GAC 548 29 Precip intensa pos-seca 2.51 Sim
ESTRUTURAS &#124; DEFICIÊNCIA EM OUTRAS INSTALAÇÕES FIXAS 258 29 Temp max >= 40°C 2.35 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA EM SINAL 21974 57 Vento >= 60 km/h 2.29 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA EM SINAL 21974 57 Chuva + vento forte 2.29 Sim
INSTALAÇÕES DE SINALIZAÇÃO EM PASSAGEM DE NÍVEL &#124; AVARIA EM PN 22372 35 Vento >= 60 km/h 2.25 Sim
INSTALAÇÕES DE SINALIZAÇÃO EM PASSAGEM DE NÍVEL &#124; AVARIA EM PN 22372 35 Chuva + vento forte 2.25 Sim
RISCOS / ACIDENTES / INCIDENTES PERIGOSOS &#124; ATOS DE VANDALISMO 7527 48 Vento >= 60 km/h 2.23 Sim
RISCOS / ACIDENTES / INCIDENTES PERIGOSOS &#124; ATOS DE VANDALISMO 7527 48 Chuva + vento forte 2.23 Sim
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; TELECOMUNICAÇÕES 4419 44 Intensidade > 10mm/h 1.95 Sim
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; SIP 3840 27 Intensidade > 10mm/h 1.73 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA EM CIRCUITO VIA 5678 50 Precip > 20mm 1.71 Sim
RISCOS / ACIDENTES / INCIDENTES PERIGOSOS &#124; FURTOS / ROUBOS / ASSALTOS 1384 47 Intensidade > 10mm/h 1.70 Sim
OUTRAS CAUSAS &#124; AVARIA EM MATERIAL CIRCULANTE 285 27 Precip > 20mm 1.68 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; AVARIA/ANOMALIA EM CONTADOR DE EIXOS 4183 31 Intensidade > 10mm/h 1.59 Sim
PESSOAL 464 28 Onda de calor 1.57 Sim
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO 511 24 Intensidade > 10mm/h 1.53 Sim
INSTALAÇÕES DE SINALIZAÇÃO &#124; CONVEL (IF) 43273 44 Onda de calor 1.45 Sim
IRREGULARIDADES NA EXECUÇÃO DE OBRAS 4736 51 Precip > 20mm 1.43 Sim
INSTALAÇÕES DE ALIMENTAÇÃO ELÉCTRICA &#124; SOBRETENSÃO 290 14 Onda de calor 1.42 Sim
INSTALAÇÕES DE TELECOMUNICAÇÕES 670 33 Precip intensa pos-seca 1.37 Sim
PEDIDO DA EMPRESA FERROVIÁRIA 4068 27 Precip intensa pos-seca 1.35 Sim
ORGANIZAÇÃO DO HORÁRIO TÉCNICO 2296 25 Precip intensa pos-seca 1.30
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; PII (PROGRAMAÇÃO INFORMÁTICA ITINERÁRIOS) 3149 24 Onda de frio 1.25
ERROS NOS PROCEDIMENTOS OPERACIONAIS 2781 43 Precip intensa pos-seca 1.24
OUTRAS CAUSAS 2398 34 Precip intensa pos-seca 1.24
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; RADIO SOLO-COMBOIO 2941 32 Onda de frio 1.18
FALHA NOS SISTEMAS DE APOIO À EXPLORAÇÃO &#124; SATA 213 17 Onda de calor 1.18
OBRAS PLANEADAS 230 26 Humidade > 90% 1.16

Relatório gerado automaticamente em 06/03/2026 22:36 · Fase 4 EDA · Plano de Trabalhos 2026

9 Exportar para PDF

# Instalar dependências se necessário
pkgs_necessarios <- c("rmarkdown", "tinytex", "knitr",
                       "sf", "rnaturalearth", "rnaturalearthdata")
pkgs_faltam <- pkgs_necessarios[!pkgs_necessarios %in% installed.packages()]
if (length(pkgs_faltam) > 0) install.packages(pkgs_faltam)

# Instalar TinyTeX se não tiver LaTeX instalado
if (!tinytex::is_tinytex() && Sys.which("xelatex") == "") {
  tinytex::install_tinytex()
  tinytex::tlmgr_install(c("booktabs", "longtable", "float",
                            "colortbl", "fancyhdr", "xcolor"))
}

rmd_path <- r"(\\refer.pt\rforfs01\DAM\AM-PR\AM-INF\Dados Climáticos Vs Ocorrências\eGOC_processed\eda_ocorrencias_clima.Rmd)"

# PDF
rmarkdown::render(
  input         = rmd_path,
  output_format = "pdf_document",
  output_file   = "eda_ocorrencias_clima.pdf",
  output_dir    = dirname(rmd_path),
  envir         = new.env()
)

# HTML (mais interactivo — tabs, TOC flutuante)
rmarkdown::render(
  input         = rmd_path,
  output_format = "html_document",
  output_file   = "eda_ocorrencias_clima.html",
  output_dir    = dirname(rmd_path),
  envir         = new.env()
)

message("Exportação concluída.")