3.1 Fonte e contexto

Os dados vêm do Kaggle (Global Earthquakes Seismic Telemetry), que reúne 10.000 eventos com magnitude entre 4,5 e 7,8. Os campos de medição (tempo, coordenadas, profundidade, magnitude, estações etc.) seguem o formato do catálogo do USGS. Colunas como magnitude_severity_class, focal_depth_regime, tsunami_hazard_potential, tectonic_boundary_zone, evidence_span, cot_rationale e confidence_score são anotações adicionadas depois à medição original.

Peculiaridades:

• time está em milissegundos desde 1970 (época Unix), não em data legível.
• place é texto livre (“45 km SW of Cidade, País”), com formatos diferentes para eventos no mar.
• As variáveis de qualidade da medição (nst, gap, dmin, rms) vêm preenchidas em todos os eventos, ou seja, não há valores ausentes para tratar.
• As classes de severidade são muito desbalanceadas (94% moderados).

eq_raw <- read_csv("~/Downloads/archive/global_earthquakes_seismic.csv")
dim(eq_raw)

## [1] 10000    25

3.2 Limpeza e transformação

A tabela abaixo mostra as primeiras linhas das colunas que mais precisam de ajuste: time (número gigante, não é uma data), place (texto livre) e as colunas de texto mag_type, net e status, que podem misturar maiúsculas e minúsculas.

# Função auxiliar para tabelas compactas
tabela <- function(df, cap) {
  df |>
    kbl(caption = cap) |>
    kable_styling(bootstrap_options = c("striped", "condensed"), full_width = FALSE, font_size = 12)
}

# Cópia só com nomes padronizados, para comparar antes/depois
eq_raw_l <- clean_names(eq_raw)

eq_raw_l |>
  dplyr::select(type, time, place, mag_type, net, status, tsunami) |>
  head(5) |>
  mutate(time = format(time, scientific = FALSE)) |>
  tabela("ANTES: primeiras linhas dos dados brutos (colunas a ajustar)")

As etapas de limpeza são:

1. Padronizar nomes das colunas (clean_names), por exemplo magType vira mag_type.
2. Filtrar apenas type == "earthquake", para não misturar outros tipos de evento.
3. Converter time de milissegundos (época Unix) para data e hora, e extrair ano, mês e hora.
4. Padronizar o texto de mag_type, status e net para minúsculas.
5. Converter para fator as classes de severidade e de profundidade, com a ordem correta dos níveis.
6. Separar place em distância, direção, local e região (regex).
7. Unir ano e mês em uma chave ano_mes.

eq <- eq_raw_l |>
  filter(type == "earthquake") |>                    # mantém só terremotos
  mutate(
    datetime = as_datetime(time / 1000),             # ms -> data/hora UTC
    ano = year(datetime), mes = month(datetime, label = TRUE), hora = hour(datetime),
    across(c(mag_type, status, net), str_to_lower),  # padroniza texto
    mag_class = factor(magnitude_severity_class,
                       levels = c("moderate", "strong", "major", "great_catastrophic")),
    depth_regime = factor(focal_depth_regime,
                          levels = c("shallow_crustal", "intermediate_subduction", "deep_mantle_slab")),
    tsunami = as.logical(tsunami)
  ) |>
  # SEPARAR: place em distância, direção, local e região
  tidyr::extract(place, into = c("dist_km", "direcao", "local", "regiao"),
          regex = "^(\\d+)\\s*km\\s+([NSEW]{1,3})\\s+of\\s+(.*?),\\s*(.*)$",
          remove = FALSE, convert = TRUE) |>
  # UNIR: chave ano-mês
  unite("ano_mes", ano, mes, sep = "-", remove = FALSE)

Veja o resultado nas mesmas linhas de antes, agora com as novas colunas geradas:

eq |>
  dplyr::select(datetime, place, dist_km, direcao, local, regiao, ano_mes, mag_type, net, status, tsunami) |>
  head(6) |>
  tabela("")

Os eventos cujo place não segue o padrão “N km DIREÇÃO of Local, Região” (por exemplo, “Mid-Atlantic Ridge”) ficam com NA nas colunas separadas, o que é esperado.

3.3 Identificando o país de cada evento

Usamos as coordenadas, que são mais confiáveis que o texto de place. Eventos no mar são atribuídos ao país mais próximo.

sf_use_s2(FALSE)  # evita problemas de geometria no join (aproximação plana)
paises <- ne_countries(scale = "medium", returnclass = "sf") |>
  dplyr::select(iso3 = adm0_a3, pais = name)

eq_sf <- st_as_sf(eq, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
eq_sf$pais_idx <- st_nearest_feature(eq_sf, paises)
eq_sf$iso3 <- paises$iso3[eq_sf$pais_idx]
eq_sf$pais <- paises$pais[eq_sf$pais_idx]

Duas colunas novas foram criadas: pais e iso3 (código de 3 letras, usado para cruzar com o Banco Mundial). Uma amostra de eventos mostra como o texto original de place se compara ao país atribuído, e a segunda tabela mostra os países com mais eventos.

set.seed(42)
eq_sf |>
  st_drop_geometry() |>
  dplyr::select(place, latitude, longitude, pais, iso3) |>
  slice_sample(n = 8) |>
  tabela("Amostra de eventos com o país atribuído pelas coordenadas")

eq_sf |>
  st_drop_geometry() |>
  count(pais, sort = TRUE, name = "eventos") |>
  head(8) |>
  tabela("Países com mais eventos")

3.4 Mesclando com dados externos

(a) Indicadores por país (Banco Mundial). Baixamos população total e PIB per capita de 2020 para cada país e juntamos pela coluna iso3.

wdi <- WDI(country = "all", indicator = c(pop_pais = "SP.POP.TOTL", pib_pc = "NY.GDP.PCAP.CD"),
           start = 2020, end = 2020) |>
  dplyr::select(iso3 = iso3c, pop_pais, pib_pc)

eq_sf <- eq_sf |> left_join(wdi, by = "iso3")

eq_sf |>
  st_drop_geometry() |>
  distinct(pais, iso3, pop_pais, pib_pc) |>
  slice_sample(n = 8) |>
  mutate(pop_pais = comma(pop_pais), pib_pc = comma(round(pib_pc))) |>
  tabela("Amostra de países após o merge com o Banco Mundial (2020)")

Após o merge, 788 eventos ficaram sem indicador do Banco Mundial (território sem código correspondente ou sem dado em 2020).

(b) População exposta ao redor do epicentro. Estimamos quantas pessoas vivem em um raio que cresce com a magnitude (15 km para M4,5, dobrando a cada unidade de magnitude). É uma heurística simples de “área sentida”, e não um modelo físico de dano.

tibble(Magnitude = c(4.5, 5.5, 6.5, 7.5)) |>
  mutate(`Raio considerado (km)` = calc_raio_km(Magnitude)) |>
  tabela("Raio de exposição por magnitude")

3.5 Conjunto final e resumo

pop <- eq_final$pop_exposta
fmt <- function(x, acc = 1) number(x, accuracy = acc, big.mark = ".", decimal.mark = ",")

tribble(
  ~Indicador, ~Valor,
  "Eventos", fmt(nrow(eq_final)),
  "Período", paste(min(as_date(eq_final$datetime)), "a", max(as_date(eq_final$datetime))),
  "Magnitude média", fmt(mean(eq_final$mag), 0.01),
  "Profundidade mediana (km)", fmt(median(eq_final$depth), 0.1),
  "Países (mais próximo do epicentro)", fmt(n_distinct(eq_final$pais)),
  "Eventos sem população no raio estimado", percent(mean(pop == 0), accuracy = 0.1, decimal.mark = ","),
  "Eventos fora da área coberta pelo raster", fmt(sum(eq_final$pop_sem_dado)),
  "População exposta: mediana (só eventos com população)", fmt(median(pop[pop > 0])),
  "População exposta: média", fmt(mean(pop)),
  "População exposta: máximo", fmt(max(pop))
) |>
  tabela("Resumo do conjunto de dados limpo")

O conjunto final tem 9998 eventos em 114 países. A mediana de população exposta é 0 porque 69,6% dos eventos ocorreram em locais sem população dentro do raio estimado, principalmente oceanos e fossas oceânicas. Considerando só os eventos que atingem alguma população, a mediana sobe para 634, e a média geral (1.739) é puxada por poucos eventos que atingem milhões de pessoas.

A tabela abaixo mostra uma amostra do conjunto final, com as principais variáveis usadas nas análises:

set.seed(2026)
eq_final |>
  dplyr::select(datetime, mag, depth, mag_class, tectonic_boundary_zone, pais, tsunami, pop_exposta) |>
  slice_sample(n = 10) |>
  mutate(datetime = format(datetime, "%Y-%m-%d %H:%M"),
         pop_exposta = comma(round(pop_exposta))) |>
  tabela("Conjunto final: amostra de 10 eventos")

4.1 Poucos eventos grandes, muitos pequenos

p <- eq_final |>
  count(mag_round = round(mag, 1)) |>
  ggplot(aes(mag_round, n)) +
  geom_col(fill = "#c0392b") +
  scale_y_log10(labels = comma) +
  labs(title = "Cada degrau de magnitude é muito mais raro",
       subtitle = "Escala logarítmica no eixo y (relação de Gutenberg-Richter)",
       x = "Magnitude", y = "Número de eventos (log)")
ggplotly(p)

4.2 Onde e a que profundidade acontecem

mundo <- ne_countries(returnclass = "sf")
ggplot() +
  geom_sf(data = mundo, fill = "grey95", color = "grey80") +
  geom_point(data = eq_final, aes(longitude, latitude, color = tectonic_boundary_zone),
             size = 0.6, alpha = 0.6) +
  labs(title = "O Anel de Fogo do Pacífico concentra a maior parte dos eventos",
       color = "Zona tectônica") +
  theme_minimal()

ggplot(eq_final, aes(depth, tectonic_boundary_zone, fill = tectonic_boundary_zone)) +
  geom_density_ridges(alpha = 0.7, show.legend = FALSE) +
  labs(title = "Eventos profundos concentram-se em zonas de subducção",
       x = "Profundidade (km)", y = NULL) +
  theme_minimal()

# FORMATO LARGO: contagem de eventos por zona e classe de magnitude
eq_final |>
  count(tectonic_boundary_zone, mag_class) |>
  pivot_wider(names_from = mag_class, values_from = n, values_fill = 0) |>
  kbl(caption = "Eventos por zona tectônica e classe de magnitude") |>
  kable_styling(full_width = FALSE)

4.3 Magnitude não é perigo: a população muda tudo

p2 <- eq_final |>
  ggplot(aes(mag, pop_exposta, color = tectonic_boundary_zone,
             text = paste0(title, "<br>Pop. exposta: ", comma(pop_exposta)))) +
  geom_point(alpha = 0.4) +
  scale_y_log10(labels = comma) +
  labs(title = "Eventos de mesma magnitude expõem populações muito diferentes",
       x = "Magnitude", y = "População exposta (log)", color = "Zona")
ggplotly(p2, tooltip = "text")

4.4 Os eventos mais fortes não são os mais perigosos

top_mag <- eq_final |> slice_max(mag, n = 10, with_ties = FALSE) |>
  dplyr::select(title, mag, pais, pop_exposta)
top_pop <- eq_final |> slice_max(pop_exposta, n = 10, with_ties = FALSE) |>
  dplyr::select(title, mag, pais, pop_exposta)

top_mag |> mutate(pop_exposta = comma(pop_exposta)) |>
  kbl(caption = "Top 10 por magnitude") |> kable_styling(full_width = FALSE)

top_pop |> mutate(pop_exposta = comma(pop_exposta)) |>
  kbl(caption = "Top 10 por população exposta") |> kable_styling(full_width = FALSE)

top_leaf <- eq_final |> slice_max(pop_exposta, n = 300, with_ties = FALSE)
pal <- colorNumeric("YlOrRd", top_leaf$mag)
leaflet(top_leaf) |>
  addTiles() |>
  addCircleMarkers(~longitude, ~latitude, radius = ~log_pop * 1.5, color = ~pal(mag),
                   fillOpacity = 0.6, stroke = FALSE,
                   popup = ~paste0(title, "<br>Mag: ", mag, "<br>Pop. exposta: ", comma(pop_exposta))) |>
  addLegend("bottomright", pal = pal, values = ~mag, title = "Magnitude")

4.5 Países: onde o risco se acumula

eq_final |>
  group_by(pais) |>
  summarise(eventos = n(), mag_max = max(mag),
            pop_exposta_total = sum(pop_exposta),
            pib_pc = first(pib_pc), .groups = "drop") |>
  arrange(desc(pop_exposta_total)) |>
  datatable(caption = "Países por população exposta acumulada", rownames = FALSE) |>
  formatRound("pop_exposta_total", digits = 0, mark = ".")

4.6 Tsunami: o que diferencia os alertas?

eq_final |>
  group_by(tsunami) |>
  summarise(n = n(), mag_mediana = median(mag), prof_mediana = median(depth), .groups = "drop") |>
  kbl(caption = "Eventos com e sem alerta de tsunami") |> kable_styling(full_width = FALSE)

# Modelo simples (opcional): o que prediz o alerta?
modelo <- glm(tsunami ~ mag + depth, data = eq_final, family = binomial)
summary(modelo)$coefficients |> round(3) |> kbl() |> kable_styling(full_width = FALSE)

4.7 Qualidade das medições (formato longo)

# FORMATO LONGO: várias métricas de qualidade em uma só coluna
eq_final |>
  dplyr::select(status, nst, gap, rms) |>
  pivot_longer(c(nst, gap, rms), names_to = "metrica", values_to = "valor") |>
  filter(!is.na(valor)) |>
  ggplot(aes(status, valor, fill = status)) +
  geom_boxplot(show.legend = FALSE) +
  facet_wrap(~ metrica, scales = "free_y") +
  labs(title = "Qualidade da localização: eventos revisados vs. automáticos",
       x = NULL, y = NULL)

4.8 O que dizem as anotações de texto

eq_final |>
  dplyr::select(mag_class, cot_rationale) |>
  unnest_tokens(palavra, cot_rationale) |>
  anti_join(stop_words, by = c("palavra" = "word")) |>
  filter(!str_detect(palavra, "^[0-9.]+$")) |>
  count(mag_class, palavra) |>
  group_by(mag_class) |>
  slice_max(n, n = 8, with_ties = FALSE) |>
  ggplot(aes(reorder_within(palavra, n, mag_class), n, fill = mag_class)) +
  geom_col(show.legend = FALSE) +
  scale_x_reordered() +
  coord_flip() +
  facet_wrap(~ mag_class, scales = "free") +
  labs(title = "Palavras mais frequentes nas justificativas, por classe", x = NULL, y = "Frequência")