Informe Descriptivo de la Sismicidad en Colombia

Análisis completo de los últimos 20 años (fuente: catálogo sísmico USGS)

# -----------------------------------------------------------------------
# Función de lectura robusta:
#   1. Intenta leer el csv de forma normal con read_csv.
#   2. Si eso falla o produce muy pocas columnas (señal típica de un
#      archivo exportado con cada fila envuelta en comillas dobles
#      adicionales), repara el archivo línea por línea y lo vuelve a leer.
# Esto permite que el mismo código funcione tanto con el csv "limpio"
# descargado directamente del USGS como con exportaciones mal formateadas.
# -----------------------------------------------------------------------
leer_datos <- function(ruta) {
  intento <- tryCatch(
    suppressWarnings(readr::read_csv(ruta, show_col_types = FALSE)),
    error = function(e) NULL
  )
  if (!is.null(intento) && ncol(intento) >= 10) {
    return(as.data.frame(intento))
  }
  raw <- readLines(ruta, encoding = "UTF-8", warn = FALSE)
  raw <- sub("^\ufeff", "", raw)                # quita BOM si existe
  encabezado <- raw[1]
  cuerpo <- raw[-1]
  cuerpo <- cuerpo[nzchar(trimws(cuerpo))]
  tiene_comillas <- startsWith(cuerpo, "\"") & endsWith(cuerpo, "\"")
  cuerpo[tiene_comillas] <- substr(cuerpo[tiene_comillas], 2,
                                    nchar(cuerpo[tiene_comillas]) - 1)
  cuerpo <- gsub("\"\"", "\"", cuerpo, fixed = TRUE)
  texto_reparado <- paste(c(encabezado, cuerpo), collapse = "\n")
  read.csv(text = texto_reparado, stringsAsFactors = FALSE, encoding = "UTF-8")
}

# Si el archivo indicado en params$archivo no existe, se busca cualquier
# .csv en la carpeta como respaldo (útil si el archivo fue renombrado).
ruta_datos <- params$archivo
if (!file.exists(ruta_datos)) {
  candidatos <- list.files(pattern = "\\.csv$", ignore.case = TRUE)
  if (length(candidatos) >= 1) ruta_datos <- candidatos[1]
}

sismos_raw <- leer_datos(ruta_datos)

# ---- Limpieza y variables derivadas ------------------------------------
sismos <- sismos_raw %>%
  mutate(
    fecha      = ymd_hms(time, tz = "UTC"),
    anio       = year(fecha),
    mes        = factor(month(fecha), levels = 1:12, labels = meses_es),
    hora       = hour(fecha),
    dia_semana = wday(fecha, label = TRUE, abbr = TRUE, week_start = 1),
    mag        = as.numeric(mag),
    depth      = as.numeric(depth),
    latitude   = as.numeric(latitude),
    longitude  = as.numeric(longitude)
  ) %>%
  filter(!is.na(fecha), !is.na(mag), !is.na(depth)) %>%
  # el catálogo incluye eventos regionales (Venezuela, Ecuador, Perú, Panamá);
  # nos quedamos únicamente con sismos localizados en Colombia
  filter(str_detect(place, "Colombia")) %>%
  # ventana de análisis: últimos N años, contados desde la fecha en que se
  # genera el informe (así el informe queda vigente si se vuelve a knittear
  # más adelante con un catálogo actualizado)
  filter(fecha >= (Sys.time() - years(params$anios_analisis))) %>%
  mutate(
    magnitud_cat = cut(
      mag,
      breaks = c(-Inf, 3, 4, 5, 6, 7, Inf),
      labels = c("Micro (<3.0)", "Menor (3.0-3.9)", "Ligero (4.0-4.9)",
                 "Moderado (5.0-5.9)", "Fuerte (6.0-6.9)", "Mayor o gran (>=7.0)"),
      right = FALSE
    ),
    profundidad_cat = cut(
      depth,
      breaks = c(-Inf, 70, 300, Inf),
      labels = c("Superficial (0-70 km)", "Intermedio (70-300 km)", "Profundo (>300 km)"),
      right = TRUE
    ),
    lugar_cercano = place %>%
      str_remove("^\\d+\\s*km\\s*[A-Za-z]+\\s*of\\s*") %>%
      str_remove(",?\\s*Colombia$") %>%
      str_trim()
  ) %>%
  arrange(fecha)

n_total       <- nrow(sismos)
n_original    <- sum(str_detect(sismos_raw$place, "Colombia"), na.rm = TRUE)
rango_fecha   <- range(sismos$fecha)
mag_min_datos <- min(sismos$mag)

0.1 Introducción

Este informe presenta una caracterización descriptiva y un complemento estadístico-inferencial de la sismicidad ocurrida en el territorio colombiano durante los últimos 20 años (06 Sep 2006 – 10 Aug 2026), con base en el catálogo sísmico del USGS (sismos_colombia.csv), filtrado a los eventos localizados en Colombia. En total se analizan 657 sismos.

Dos precisiones importantes sobre los datos, necesarias para interpretar correctamente los resultados:

1. La consulta original al USGS trae un radio alrededor de Colombia (por eso el archivo también contiene sismos de Venezuela, Ecuador, Perú y Panamá); de los 1285 registros del archivo, 784 están ubicados en Colombia, y de esos, 657 caen dentro de la ventana de 20 años analizada.
2. La magnitud mínima registrada en el archivo es 4.5, es decir, la base no incluye sismos menores (magnitud < 4.5). Por eso, en las tablas de clasificación por magnitud aparecerán categorías como “Micro” o “Menor” con 0 (o muy pocos) eventos: no es que esos sismos no ocurran en Colombia, sino que este catálogo específico no los reporta.

El documento se organiza en cuatro capítulos: (1) descripción por rangos de magnitud y profundidad, (2) descripción temporal, (3) descripción espacial y (4) conclusiones y recomendaciones. En cada capítulo, además de las tablas y gráficos descriptivos, se incluye un análisis estadístico complementario (pruebas de hipótesis, regresión, agrupamiento) que permite sustentar con mayor rigor las conclusiones.

---

1 Capítulo 1. Descripción de los sismos por rangos de magnitud y profundidad

1.1 Estadísticas generales

skew <- function(x) {
  n <- length(x); m <- mean(x); s <- sd(x)
  (sum((x - m)^3) / n) / s^3
}
kurt <- function(x) {
  n <- length(x); m <- mean(x); s <- sd(x)
  (sum((x - m)^4) / n) / s^4 - 3
}

sismos %>%
  summarise(
    `N° sismos`               = n(),
    `Magnitud mínima`         = round(min(mag), 1),
    `Magnitud media`          = round(mean(mag), 2),
    `Magnitud mediana`        = round(median(mag), 2),
    `Desv. estándar magnitud` = round(sd(mag), 2),
    `Asimetría magnitud`      = round(skew(mag), 2),
    `Curtosis magnitud`       = round(kurt(mag), 2),
    `Magnitud máxima`         = round(max(mag), 1),
    `Profundidad media (km)`  = round(mean(depth), 1),
    `Profundidad mediana (km)` = round(median(depth), 1),
    `Profundidad máxima (km)` = round(max(depth), 1)
  ) %>%
  pivot_longer(everything(), names_to = "Estadístico", values_to = "Valor") %>%
  kable(caption = "Estadísticas descriptivas del período analizado") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Estadísticas descriptivas del período analizado

Estadístico	Valor
N° sismos	657.00
Magnitud mínima	4.50
Magnitud media	4.75
Magnitud mediana	4.60
Desv. estándar magnitud	0.36
Asimetría magnitud	3.09
Curtosis magnitud	13.87
Magnitud máxima	7.40
Profundidad media (km)	85.40
Profundidad mediana (km)	64.40
Profundidad máxima (km)	436.50

La asimetría (skewness) de la magnitud es positiva, lo que confirma numéricamente algo esperado en sismología: hay muchos más sismos de magnitud baja que de magnitud alta, con una “cola” larga hacia la derecha (pocos eventos muy grandes).

1.2 Distribución por rango de magnitud

tabla_mag <- sismos %>%
  count(magnitud_cat, .drop = FALSE, name = "n_sismos") %>%
  mutate(porcentaje = percent(n_sismos / sum(n_sismos), accuracy = 0.1))

tabla_mag %>%
  kable(col.names = c("Rango de magnitud", "N° de sismos", "% del total"),
        caption = "Sismos por rango de magnitud") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Sismos por rango de magnitud

Rango de magnitud	N° de sismos	% del total
Micro (<3.0)	0	0.0%
Menor (3.0-3.9)	0	0.0%
Ligero (4.0-4.9)	538	81.9%
Moderado (5.0-5.9)	106	16.1%
Fuerte (6.0-6.9)	11	1.7%
Mayor o gran (>=7.0)	2	0.3%

ggplot(tabla_mag, aes(x = magnitud_cat, y = n_sismos, fill = magnitud_cat)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = n_sismos), vjust = -0.4, size = 3.5) +
  scale_fill_manual(values = paleta_mag, drop = FALSE) +
  scale_x_discrete(drop = FALSE) +
  labs(title = "Número de sismos por rango de magnitud",
       x = "Rango de magnitud", y = "Número de sismos") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

ggplot(sismos, aes(x = mag)) +
  geom_histogram(binwidth = 0.2, boundary = 4.5, fill = "#e67e22", color = "white") +
  labs(title = "Distribución continua de la magnitud",
       x = "Magnitud", y = "Número de sismos") +
  theme_minimal(base_size = 12)

1.3 Distribución por rango de profundidad

tabla_prof <- sismos %>% count(profundidad_cat, .drop = FALSE, name = "n_sismos")

ggplot(tabla_prof, aes(x = profundidad_cat, y = n_sismos, fill = profundidad_cat)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = n_sismos), vjust = -0.4, size = 3.5) +
  scale_fill_manual(values = c("#3498db", "#f39c12", "#c0392b")) +
  labs(title = "Número de sismos por rango de profundidad",
       x = "Rango de profundidad", y = "Número de sismos") +
  theme_minimal(base_size = 12)

ggplot(sismos, aes(x = profundidad_cat, y = depth, fill = profundidad_cat)) +
  geom_boxplot(show.legend = FALSE) +
  scale_fill_manual(values = c("#3498db", "#f39c12", "#c0392b")) +
  labs(title = "Profundidad de los sismos, según rango",
       x = "Rango de profundidad", y = "Profundidad (km)") +
  theme_minimal(base_size = 12)

1.4 Relación magnitud – profundidad

sismos %>%
  count(magnitud_cat, profundidad_cat, .drop = FALSE) %>%
  ggplot(aes(x = magnitud_cat, y = profundidad_cat, fill = n)) +
  geom_tile(color = "white") +
  geom_text(aes(label = n), color = "white", size = 3.5) +
  scale_fill_gradient(low = "#a1d99b", high = "#08519c") +
  scale_x_discrete(drop = FALSE) +
  labs(title = "Cruce entre rango de magnitud y rango de profundidad",
       x = "Rango de magnitud", y = "Rango de profundidad", fill = "N° sismos") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

ggplot(sismos, aes(x = depth, y = mag, color = magnitud_cat)) +
  geom_point(alpha = 0.6) +
  scale_color_manual(values = paleta_mag, drop = FALSE) +
  labs(title = "Dispersión magnitud vs. profundidad",
       x = "Profundidad (km)", y = "Magnitud", color = "Rango de magnitud") +
  theme_minimal(base_size = 12)

1.4.1 Complemento estadístico: ¿magnitud y profundidad están asociadas?

test_cor <- cor.test(sismos$mag, sismos$depth, method = "pearson")

Se calcula el coeficiente de correlación de Pearson entre magnitud y profundidad, junto con una prueba de hipótesis (\(H_0\): la correlación poblacional es 0):

• Coeficiente de correlación: r = -0.098
• Valor p: 0.0121

Con un nivel de significancia del 5%, la correlación es estadísticamente significativa, aunque su magnitud (r = -0.1) indica una asociación lineal débil entre magnitud y profundidad: la profundidad por sí sola no permite predecir bien la magnitud de un sismo.

Adicionalmente, se evalúa si la magnitud y el rango de profundidad son independientes mediante una prueba chi-cuadrado sobre la tabla de contingencia mostrada arriba:

tabla_contingencia <- table(sismos$magnitud_cat, sismos$profundidad_cat)
# se eliminan categorías de magnitud sin observaciones (estructuralmente
# vacías por el límite inferior de la base, ver Introducción) para que la
# prueba chi-cuadrado sea válida
tabla_contingencia <- tabla_contingencia[rowSums(tabla_contingencia) > 0, , drop = FALSE]
test_chi2 <- suppressWarnings(chisq.test(tabla_contingencia))

• Estadístico chi-cuadrado: 21.34
• Grados de libertad: 6
• Valor p: 0.00159

El valor p es menor a 0.05, por lo que se rechaza la hipótesis de independencia: el rango de profundidad sí está asociado con el rango de magnitud del sismo.

Nota: dado que algunas celdas de la tabla tienen pocas observaciones, este resultado debe interpretarse con cautela (la prueba chi-cuadrado es más confiable cuando las frecuencias esperadas son mayores o iguales a 5).

1.5 Los 10 sismos de mayor magnitud del período

sismos %>%
  slice_max(mag, n = 10) %>%
  transmute(
    Fecha = format(fecha, "%Y-%m-%d"),
    Lugar = lugar_cercano,
    Magnitud = mag,
    `Profundidad (km)` = round(depth, 1)
  ) %>%
  kable(caption = "Diez sismos de mayor magnitud registrados") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Diez sismos de mayor magnitud registrados

Fecha	Lugar	Magnitud	Profundidad (km)
2026-08-10	San José del Palmar	7.4	110.3
2012-09-30	San Agustín	7.3	170.0
2013-02-09	Yacuanquer	6.9	145.0
2007-09-10	Timbiquí	6.8	15.0
2013-08-13	Bahía Solano	6.7	12.0
2025-06-08	Paratebueno	6.3	9.0
2007-03-18	Santa Genoveva de Docordó	6.2	7.0
2015-03-10	Cepitá	6.2	155.0
2019-03-23	Versalles	6.1	122.0
2023-08-17	El Calvario	6.1	10.0

1.6 Complemento: ley de Gutenberg-Richter (relación frecuencia-magnitud)

Una de las relaciones más conocidas en sismología es la ley de Gutenberg-Richter, que dice que el logaritmo del número de sismos con magnitud igual o mayor a M decrece linealmente con M:

\[\log_{10} N(\geq M) = a - bM\]

El parámetro b (pendiente cambiada de signo) se estima aquí mediante una regresión lineal simple (mínimos cuadrados), un método visto en el curso, aplicado sobre el logaritmo de la frecuencia acumulada:

pasos_mag <- seq(floor(min(sismos$mag) * 10) / 10,
                  ceiling(max(sismos$mag) * 10) / 10 - 0.5, by = 0.1)

frec_acum <- tibble(magnitud = pasos_mag) %>%
  rowwise() %>%
  mutate(n_acum = sum(sismos$mag >= magnitud)) %>%
  ungroup() %>%
  filter(n_acum >= 5) %>%          # se excluye la cola con muy pocos datos
  mutate(log_n = log10(n_acum))

modelo_gr <- lm(log_n ~ magnitud, data = frec_acum)
b_value   <- round(-coef(modelo_gr)[["magnitud"]], 2)
a_value   <- round(coef(modelo_gr)[["(Intercept)"]], 2)
r2_gr     <- round(summary(modelo_gr)$r.squared, 3)

ggplot(frec_acum, aes(x = magnitud, y = log_n)) +
  geom_point(color = "#2c3e50") +
  geom_smooth(method = "lm", se = FALSE, color = "firebrick") +
  labs(title = "Relación frecuencia-magnitud (ley de Gutenberg-Richter)",
       subtitle = paste0("b = ", b_value, "  |  a = ", a_value, "  |  R² = ", r2_gr),
       x = "Magnitud (M)", y = expression(log[10]~"N(">=M*")")) +
  theme_minimal(base_size = 12)

El valor estimado es b = 0.98 (con R² = 0.965, un ajuste muy bueno). Valores de b cercanos a 1 son típicos a nivel global; un valor cercano a 1, en línea con el comportamiento típico de la sismicidad regional. Este parámetro suele usarse como insumo para estudios de amenaza sísmica.

---

2 Capítulo 2. Descripción temporal de los sismos

2.1 Frecuencia de ocurrencia según la intensidad

La siguiente tabla estima, para cada rango de magnitud, cada cuántos días ocurre en promedio un sismo de esa intensidad en Colombia (a partir del tiempo transcurrido entre eventos consecutivos de la misma categoría).

frecuencia_cat <- sismos %>%
  arrange(magnitud_cat, fecha) %>%
  group_by(magnitud_cat, .drop = FALSE) %>%
  mutate(dias_desde_anterior = as.numeric(difftime(fecha, lag(fecha), units = "days"))) %>%
  summarise(
    n_eventos = n(),
    promedio_dias_entre_eventos = round(mean(dias_desde_anterior, na.rm = TRUE), 1),
    mediana_dias_entre_eventos  = round(median(dias_desde_anterior, na.rm = TRUE), 1),
    .groups = "drop"
  )

frecuencia_cat %>%
  kable(col.names = c("Rango de magnitud", "N° eventos", "Días prom. entre eventos",
                       "Mediana días entre eventos"),
        caption = "Frecuencia de ocurrencia por rango de magnitud") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Frecuencia de ocurrencia por rango de magnitud

Rango de magnitud	N° eventos	Días prom. entre eventos	Mediana días entre eventos
Micro (<3.0)	0	NaN	NA
Menor (3.0-3.9)	0	NaN	NA
Ligero (4.0-4.9)	538	13.5	8.9
Moderado (5.0-5.9)	106	67.9	46.1
Fuerte (6.0-6.9)	11	665.8	563.7
Mayor o gran (>=7.0)	2	5061.8	5061.8

frecuencia_cat %>%
  filter(!is.na(promedio_dias_entre_eventos)) %>%
  ggplot(aes(x = magnitud_cat, y = promedio_dias_entre_eventos, fill = magnitud_cat)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = promedio_dias_entre_eventos), vjust = -0.4, size = 3.5) +
  scale_fill_manual(values = paleta_mag, drop = FALSE) +
  labs(title = "Días promedio entre sismos, según rango de magnitud",
       x = "Rango de magnitud", y = "Días promedio entre eventos") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

sismos %>%
  arrange(fecha) %>%
  group_by(magnitud_cat, .drop = FALSE) %>%
  mutate(dias_desde_anterior = as.numeric(difftime(fecha, lag(fecha), units = "days"))) %>%
  ungroup() %>%
  filter(!is.na(dias_desde_anterior), dias_desde_anterior > 0) %>%
  ggplot(aes(x = magnitud_cat, y = dias_desde_anterior, fill = magnitud_cat)) +
  geom_boxplot(show.legend = FALSE) +
  scale_fill_manual(values = paleta_mag, drop = FALSE) +
  scale_y_log10() +
  labs(title = "Días entre sismos consecutivos, según rango de magnitud",
       subtitle = "Eje Y en escala logarítmica por la dispersión de los datos",
       x = "Rango de magnitud", y = "Días entre sismos (log10)") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

2.2 Serie de tiempo anual

serie_anual <- sismos %>% count(anio, name = "n_sismos")

ggplot(serie_anual, aes(x = anio, y = n_sismos)) +
  geom_col(fill = "#2c7fb8") +
  geom_smooth(method = "loess", se = FALSE, color = "firebrick", linewidth = 0.9) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(title = "Número de sismos por año en Colombia",
       subtitle = paste0("Periodo ", min(serie_anual$anio), " - ", max(serie_anual$anio)),
       x = "Año", y = "Número de sismos") +
  theme_minimal(base_size = 12)

2.2.1 Complemento estadístico: ¿hay una tendencia significativa?

modelo_tendencia <- lm(n_sismos ~ anio, data = serie_anual)
resumen_tend <- summary(modelo_tendencia)
pendiente    <- round(coef(modelo_tendencia)[["anio"]], 2)
p_tendencia  <- signif(resumen_tend$coefficients["anio", "Pr(>|t|)"], 3)

Se ajusta una regresión lineal simple (número de sismos por año en función del año) para evaluar si existe una tendencia significativa:

• Pendiente estimada: 0.91 sismos/año
• Valor p de la pendiente: 0.0224

direccion_tend <- ifelse(pendiente > 0, "creciente", "decreciente")

El valor p es menor a 0.05, por lo que la tendencia creciente es estadísticamente significativa. Sin embargo, esto puede reflejar tanto cambios reales en la actividad sísmica como mejoras en la capacidad de detección de las redes sismológicas a lo largo del período analizado, por lo que la interpretación debe ser cautelosa.

2.3 Serie anual por rango de magnitud

sismos %>%
  count(anio, magnitud_cat, .drop = FALSE) %>%
  ggplot(aes(x = anio, y = n, fill = magnitud_cat)) +
  geom_col() +
  scale_fill_manual(values = paleta_mag, drop = FALSE) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  labs(title = "Composición anual de sismos por rango de magnitud",
       x = "Año", y = "Número de sismos", fill = "Rango de magnitud") +
  theme_minimal(base_size = 12)

2.4 Estacionalidad: mes, día de la semana y hora del día

sismos %>%
  count(mes) %>%
  ggplot(aes(x = mes, y = n)) +
  geom_col(fill = "#41b6c4") +
  labs(title = "Número de sismos por mes (acumulado en todos los años)",
       x = "Mes", y = "Número de sismos") +
  theme_minimal(base_size = 12)

p_semana <- sismos %>%
  count(dia_semana) %>%
  ggplot(aes(x = dia_semana, y = n)) +
  geom_col(fill = "#88419d") +
  labs(title = "Sismos por día de la semana", x = NULL, y = "N° de sismos") +
  theme_minimal(base_size = 11)

p_hora <- ggplot(sismos, aes(x = hora)) +
  geom_histogram(binwidth = 1, fill = "#df65b0", color = "white") +
  labs(title = "Sismos por hora del día (UTC)", x = "Hora (UTC)", y = "N° de sismos") +
  theme_minimal(base_size = 11)

gridExtra::grid.arrange(p_semana, p_hora, ncol = 2)

Como es esperable desde el punto de vista físico, no se observan patrones fuertes por día de la semana ni por hora del día: la ocurrencia de sismos no depende de convenciones humanas como el calendario laboral.

2.5 Complemento: ¿los sismos se comportan como un proceso de Poisson?

En probabilidad, un supuesto común (y útil de contrastar) es que los eventos raros ocurren según un proceso de Poisson, lo que implica que el tiempo entre eventos consecutivos sigue una distribución exponencial. Se evalúa este supuesto para el conjunto completo de sismos (todas las magnitudes juntas) mediante una prueba de bondad de ajuste de Kolmogorov-Smirnov:

tiempos_entre_eventos <- sismos %>%
  arrange(fecha) %>%
  mutate(dias = as.numeric(difftime(fecha, lag(fecha), units = "days"))) %>%
  filter(!is.na(dias), dias > 0) %>%
  pull(dias)

tasa_hat <- 1 / mean(tiempos_entre_eventos)
ks_result <- suppressWarnings(
  ks.test(tiempos_entre_eventos, "pexp", rate = tasa_hat)
)

ggplot(data.frame(dias = tiempos_entre_eventos), aes(x = dias)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "#74a9cf", color = "white") +
  stat_function(fun = dexp, args = list(rate = tasa_hat),
                color = "firebrick", linewidth = 1) +
  labs(title = "Tiempo entre sismos consecutivos vs. distribución exponencial ajustada",
       subtitle = paste0("Tasa estimada = ", round(tasa_hat, 3), " sismos/día  |  ",
                          "media = ", round(mean(tiempos_entre_eventos), 1), " días"),
       x = "Días entre sismos consecutivos", y = "Densidad") +
  theme_minimal(base_size = 12)

• Estadístico D de Kolmogorov-Smirnov: 0.063
• Valor p: 0.0103

El valor p es menor a 0.05: se rechaza la hipótesis de que los tiempos entre sismos sigan una distribución exponencial pura. Esto es consistente con lo que suele observarse en catálogos sísmicos reales: los sismos no son totalmente independientes entre sí (hay réplicas y agrupamiento temporal/espacial de la actividad), por lo que un proceso de Poisson homogéneo es solo una aproximación simplificada.

---

3 Capítulo 3. Descripción espacial de los sismos

3.1 Zonas (localidades) con mayor frecuencia de sismos

top_lugares <- sismos %>%
  count(lugar_cercano, sort = TRUE) %>%
  slice_head(n = 15)

top_lugares %>%
  kable(col.names = c("Localidad de referencia", "N° de sismos"),
        caption = "15 localidades con mayor número de sismos registrados cerca") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

15 localidades con mayor número de sismos registrados cerca

Localidad de referencia	N° de sismos
Cepitá	116
Jordán	48
Nuquí	28
Aratoca	20
Bahía Solano	19
Pizarro	16
Riosucio	16
Santa Genoveva de Docordó	14
Cubará	10
Piedecuesta	10
Murindó	9
Mosquera	8
northern	8
Lejanías	7
Mesetas	7

ggplot(top_lugares, aes(x = reorder(lugar_cercano, n), y = n)) +
  geom_col(fill = "#e34a33") +
  coord_flip() +
  labs(title = "Localidades con mayor frecuencia de sismos",
       x = NULL, y = "Número de sismos") +
  theme_minimal(base_size = 12)

3.2 Agrupación de zonas de mayor actividad (clúster geográfico)

Para identificar zonas geográficas de mayor concentración sísmica más allá del nombre puntual de la localidad, se agrupan los eventos mediante un análisis de conglomerados (k-means) sobre sus coordenadas.

set.seed(123)
k_zonas <- 6

coords <- sismos %>% select(longitude, latitude) %>% drop_na()
clust  <- kmeans(coords, centers = k_zonas, nstart = 25)

sismos_geo <- sismos %>%
  drop_na(longitude, latitude) %>%
  mutate(zona_cluster = factor(clust$cluster))

resumen_zonas <- sismos_geo %>%
  group_by(zona_cluster) %>%
  summarise(
    n_sismos    = n(),
    lat_centro  = round(mean(latitude), 2),
    lon_centro  = round(mean(longitude), 2),
    mag_media   = round(mean(mag), 2),
    prof_media  = round(mean(depth), 1),
    .groups = "drop"
  ) %>%
  arrange(desc(n_sismos))

resumen_zonas %>%
  kable(col.names = c("Zona", "N° sismos", "Latitud centro", "Longitud centro",
                       "Magnitud media", "Profundidad media (km)"),
        caption = "Resumen de zonas de agrupación sísmica") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Resumen de zonas de agrupación sísmica

Zona	N° sismos	Latitud centro	Longitud centro	Magnitud media	Profundidad media (km)
2	264	6.80	-73.01	4.69	141.4
5	151	5.12	-77.14	4.81	55.4
1	81	7.44	-76.41	4.73	32.8
6	66	4.14	-73.38	4.85	27.2
4	62	2.84	-76.80	4.87	71.6
3	33	11.14	-73.48	4.70	47.0

3.2.1 Complemento estadístico: ¿la magnitud promedio difiere entre zonas?

Se aplica un análisis de varianza (ANOVA) de un factor para evaluar si la magnitud promedio es significativamente distinta entre las 6 zonas geográficas identificadas:

anova_zonas <- aov(mag ~ zona_cluster, data = sismos_geo)
resumen_anova <- summary(anova_zonas)
p_anova <- signif(resumen_anova[[1]][["Pr(>F)"]][1], 3)
f_anova <- round(resumen_anova[[1]][["F value"]][1], 2)

• Estadístico F: 5.29
• Valor p: 8.97^{-5}

El valor p es menor a 0.05: existe evidencia de que la magnitud promedio sí difiere significativamente entre al menos una de las zonas geográficas identificadas.

ggplot(sismos_geo, aes(x = zona_cluster, y = mag, fill = zona_cluster)) +
  geom_boxplot(show.legend = FALSE) +
  labs(title = "Magnitud de los sismos por zona geográfica",
       x = "Zona", y = "Magnitud") +
  theme_minimal(base_size = 12)

3.3 Mapa de distribución espacial

# Se intenta primero un mapa con fronteras administrativas obtenidas de
# rnaturalearth (más detallado, pero requiere paquetes adicionales e
# internet). Si no está disponible o falla la descarga, se usa como
# respaldo el mapa base del paquete 'maps' (viene con sus propios datos,
# por lo que siempre funciona, incluso sin conexión a internet).
colombia_sf <- NULL
if (mapa_avanzado_disponible) {
  colombia_sf <- tryCatch(
    rnaturalearth::ne_countries(country = "Colombia", scale = "medium", returnclass = "sf"),
    error = function(e) NULL
  )
}

if (!is.null(colombia_sf)) {
  ggplot() +
    geom_sf(data = colombia_sf, fill = "grey96", color = "grey40") +
    geom_point(data = sismos_geo,
               aes(x = longitude, y = latitude, color = zona_cluster, size = mag),
               alpha = 0.7) +
    coord_sf(xlim = c(-82, -66), ylim = c(-5, 13), expand = FALSE) +
    labs(title = "Distribución espacial de los sismos en Colombia",
         subtitle = "Color = zona de agrupación | Tamaño = magnitud",
         x = NULL, y = NULL, color = "Zona", size = "Magnitud") +
    theme_minimal(base_size = 12)
} else {
  ggplot(sismos_geo) +
    borders("world", regions = "Colombia", fill = "grey96", colour = "grey50") +
    geom_point(aes(x = longitude, y = latitude, color = zona_cluster, size = mag),
               alpha = 0.7) +
    coord_quickmap(xlim = c(-82, -66), ylim = c(-5, 13)) +
    labs(title = "Distribución espacial de los sismos en Colombia",
         subtitle = "Color = zona de agrupación | Tamaño = magnitud",
         x = "Longitud", y = "Latitud", color = "Zona", size = "Magnitud") +
    theme_minimal(base_size = 12)
}

3.4 Mapa de densidad (zonas “calientes”)

ggplot(sismos_geo, aes(x = longitude, y = latitude)) +
  borders("world", regions = "Colombia", fill = "grey92", colour = "grey60") +
  stat_density_2d(aes(fill = after_stat(level)), geom = "polygon", alpha = 0.5) +
  geom_point(size = 0.4, alpha = 0.4, color = "black") +
  scale_fill_viridis_c(option = "inferno") +
  coord_quickmap(xlim = c(-82, -66), ylim = c(-5, 13)) +
  labs(title = "Densidad espacial de los sismos en Colombia",
       x = "Longitud", y = "Latitud", fill = "Densidad") +
  theme_minimal(base_size = 12)

# Mapa interactivo (solo se genera si el paquete 'leaflet' está instalado).
# Al abrir el HTML sin conexión a internet, los mapas de fondo (teselas) no
# cargarán, pero los puntos y la información emergente sí funcionan.
pal <- colorFactor(palette = "Set1", domain = sismos_geo$zona_cluster)

leaflet(sismos_geo) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addCircleMarkers(
    lng = ~longitude, lat = ~latitude,
    radius = ~pmax(mag * 1.6, 3),
    color = ~pal(zona_cluster),
    stroke = FALSE, fillOpacity = 0.7,
    popup = ~paste0("<b>", lugar_cercano, "</b><br>",
                     "Fecha: ", format(fecha, "%Y-%m-%d"), "<br>",
                     "Magnitud: ", mag, "<br>",
                     "Profundidad: ", round(depth, 1), " km")
  ) %>%
  addLegend("bottomright", pal = pal, values = ~zona_cluster, title = "Zona")

---

4 Capítulo 4. Conclusiones generales y recomendaciones

zona_top      <- resumen_zonas$zona_cluster[1]
n_zona_top    <- resumen_zonas$n_sismos[1]
cat_mas_freq  <- tabla_mag$magnitud_cat[which.max(tabla_mag$n_sismos)]
prof_mas_freq <- tabla_prof$profundidad_cat[which.max(tabla_prof$n_sismos)]
anio_max      <- serie_anual$anio[which.max(serie_anual$n_sismos)]
n_anio_max    <- max(serie_anual$n_sismos)

4.1 Conclusiones generales

Descriptivas

• Durante los últimos 20 años se registraron 657 sismos de magnitud 4.5 o mayor en territorio colombiano, con magnitudes entre 4.5 y 7.4.
• La mayor parte de los eventos corresponde al rango de magnitud Ligero (4.0-4.9), consistente con la relación de Gutenberg-Richter: hay muchos más sismos moderados que sismos grandes.
• En cuanto a profundidad, predominan los sismos de categoría Superficial (0-70 km).
• La zona identificada como 2 concentra el mayor número de eventos (264 sismos), y la localidad puntual con más sismos reportados es Cepitá.
• El año 2023 registró el mayor número de sismos (54) dentro del período estudiado.
• Los sismos moderados ocurren con una frecuencia de días a pocas semanas, mientras que los sismos fuertes o mayores son mucho menos frecuentes pero de mayor impacto potencial (capítulo 2).

Estadísticas / inferenciales

• El ajuste a la ley de Gutenberg-Richter dio un valor b = 0.98 (R² = 0.965), en línea con el comportamiento típico de la sismicidad regional.
• La correlación entre magnitud y profundidad fue de r = -0.1 (p = 0.012), es decir, estadísticamente significativa pero débil.
• La prueba chi-cuadrado indica que el rango de profundidad y el rango de magnitud SÍ están asociados (p = 0.0016).
• La tendencia lineal anual resultó estadísticamente significativa (p = 0.0224), por lo que sí se observa un cambio sistemático en el número de sismos reportados por año (aunque podría deberse en parte a mejoras en el monitoreo).
• La prueba de Kolmogorov-Smirnov rechazó el supuesto de que los tiempos entre sismos siguen un proceso de Poisson puro (p = 0.01), lo que sugiere cierto agrupamiento temporal de la actividad sísmica (p. ej. réplicas).
• El ANOVA detectó diferencias significativas en la magnitud promedio entre las zonas geográficas identificadas (p = 8.97^{-5}).

4.2 Recomendaciones para la toma de decisiones

1. Priorizar el monitoreo y mantenimiento de la red sismológica en las zonas identificadas con mayor concentración de eventos (2, Cepitá), dado que son las de mayor probabilidad de actividad futura.
2. Reforzar la normativa de sismo-resistencia y la actualización de microzonificación sísmica en los municipios ubicados dentro de las zonas de mayor frecuencia, en coordinación con las autoridades locales de gestión del riesgo.
3. Mantener y divulgar planes de preparación y respuesta ante emergencias (simulacros, rutas de evacuación, protocolos escolares/hospitalarios), dado que, si bien los sismos de gran magnitud son poco frecuentes, su impacto potencial es alto.
4. Incorporar la profundidad en la evaluación de amenaza: los sismos superficiales, aun de magnitud moderada, suelen generar mayor intensidad percibida y daño que los sismos profundos de magnitud similar.
5. Usar el valor de b estimado (0.98) como insumo de referencia para estudios de amenaza sísmica probabilística, complementándolo con catálogos más completos que incluyan sismos de menor magnitud.
6. No asumir independencia total entre eventos al planear la respuesta ante un sismo importante: el resultado de la prueba de Poisson sugiere vigilar la posibilidad de réplicas o actividad agrupada en los días siguientes a un evento significativo.
7. Actualizar este informe periódicamente (por ejemplo, de forma anual) con datos recientes del catálogo sísmico, para hacer seguimiento a tendencias y detectar cambios en los patrones espacio-temporales de la sismicidad.
8. Complementar este análisis con estudios de amenaza probabilística y con la sobreposición de las zonas de mayor frecuencia sobre mapas de exposición (población, infraestructura crítica) para priorizar inversiones en reducción de riesgo.

Nota metodológica: este es un análisis descriptivo y estadístico exploratorio basado en el catálogo sísmico entregado (fuente: USGS, solo sismos M ≥ 4.5). No constituye un estudio de amenaza o riesgo sísmico formal; para decisiones normativas o de inversión se recomienda contrastar estos resultados con los estudios oficiales del Servicio Geológico Colombiano (SGC).