En este breve tutorial se muestra como cargar en R la base de datos del catálogo de sismos del Sistema Sisimológico Nacional (SSN) adscrito al Instituto de Geofísica de la UNAM. Así, como realizar algunos análisis básicos.

Servicio Sísmico Nacional

El Servicio Sísmico Nacional (SSN) de México es la institución encargada de monitorear y estudiar los sismos que ocurren en el territorio mexicano. Algunos datos y características importantes del SSN son los siguientes:

  1. Monitoreo de sismos: El SSN registra y monitorea los eventos sísmicos que ocurren en México utilizando una red de estaciones sismológicas distribuidas en todo el país.

  2. Base de datos sísmicos (catálogo de sismos): El SSN mantiene una base de datos que recopila información detallada sobre los sismos registrados, incluyendo su ubicación, magnitud, profundidad y otros parámetros asociados.

  3. Difusión de información: El SSN proporciona información en tiempo real sobre los sismos registrados a través de su sitio web y otras plataformas. Esto incluye datos como la ubicación y magnitud preliminar del sismo, así como información actualizada a medida que se obtiene.

Descargando los datos

Para descargar la base de datos vamos a la página del SSN:

http://www2.ssn.unam.mx:8080/catalogo/

Primero, se selecciona el periodo. Aquí descargamos la mayor cantidad de información posible, i.e. 1900-01-01 al último día disponible. A continuación, Magnitud (todas), Profundidad (todas), Filtrar por estado (todos). Para el formato de salida elegimos archivo de valores separados por comas. Comprobamos que no somos un robot y aceptamos las condiciones que aparecen.

Se se recomienda abrir con Excel (o con ningún programa similar a Excel), la base que se acaba de descargar. Si se quiere ver la base de datos en una hoja de cálculo, se crea una copia y abrimos la copia. El detalle es que la base de datos tiene fechas y está en csv. Excel (o programas similares) suelen modificar esta información al abrir el archivo y la función fread no será capaz de identificar este tipo de variable correctamente. Aún así seróa posible cambiar el tipo de variable, pero con esta precaución se facilitan las cosas.

Al abrir la copia de la base vemos que además de la tabla de datos, en el encabezado y al pié hay información. En el encabezado viene el doi (identificador único y permanente para las publicaciones electrónicas), se especifica que la base es del SSN, un aviso, el periodo en el que se solicitó la información y e número de eventos. Mientras que al pié de la base viene lo que sería propiamente el diccionario. Quizá sería mejor incluír esto en un diccionario de datos a parte, tal como lo hace INEGI.

Los datos

La base de datos del SSN se estructura en forma de tabla, en donde cada fila representa un evento sísmico registrado y cada columna contiene información específica sobre ese evento. Las columnas (variables) en la base de datos incluyen:

  1. Fecha y hora: Indica la fecha y hora exacta del evento sísmico.
  2. Ubicación: Proporciona la ubicación geográfica del epicentro del sismo, utilizando coordenadas de latitud y longitud.
  3. Magnitud: Indica la magnitud del sismo, que es una medida de la energía liberada durante el evento. Existe un breve documento del SSN que explica con mayor detalle la magnitud de un sismo

http://www.ssn.unam.mx/jsp/reportesEspeciales/Magnitud-de-un-sismo.pdf

aparece un warning, en particular no he tenido problemas y se descarga un documento en pdf.

  1. Profundidad: Especifica la profundidad a la que ocurrió el evento sísmico medida en kilómetros.
  2. Referencia: una referencia a una localidad importante en cuanto a numero de habitantes y cercana al epicentro.

Leyendo la historia de como fue evolucionando la medición de los sismos en México (en la pestaña ACERCA DE > HISTORIA), así como la manera en que se fue constituyendo el SSN. Es claro que en los datos observaremos gran heterogeneidad en el tiempo. En palabras simples, ahora se pueden medir mucho más sismos que hace 20, 50 o 100 años. Se detectan más sismos, se mide mejor su magnitud y ubicación. Pero eso no quiere decir que hoy en día haya más sismos que hace 100 años. Lo que nos indica es que hace 100 años no existía la capacidad técnica y tecnológica que se tiene hoy en día. Además, en los datos podemos ver ese aumento gradual en la capacidad de medición. Es importante tener en cuenta esto al analizar la información de los sismos.

Lectura de la base de datos

Para leer la base de datos usaremos la función fread de la librería data.table de R. Esta función automáticamente detacta la tabla y lee la información arrojando un warning. Esto debido a que se incluye información adicional además de la tabla. Ya que la función fread arroja una estructura de datos particular, de inmediato la convertimos a una data.frame que es fácil de manejar.

library(data.table)
my_dir <- "/Users/carloserwin/Library/CloudStorage/GoogleDrive-carloserwin@sigma.iimas.unam.mx/Mi\ unidad/IIMAS/PAGINA_WEB/tutoriales/"
DATOS0 <- data.frame(fread(paste0(my_dir, "SISMOS/SSNMX_catalogo_19000101_20230713.csv")))
## Warning in fread(paste0(my_dir, "SISMOS/SSNMX_catalogo_19000101_20230713.csv")):
## Stopped early on line 283765. Expected 10 fields but found 1. Consider fill=TRUE
## and comment.char=. First discarded non-empty line: <<"Fecha y hora local en
## tiempo del centro de Mexico. Coordenadas geograficas (latitud y longitud) del
## epicentro en grados decimales. Profundidad en kilometros.">>
nrow(DATOS0)
## [1] 283759

Se tienen aproximadamente 300 mil eventos sísmicos registrados, sin embargo la variable Magnitud se lee como texto pues los valores perdidos el SSN los codifica como no calculable (podría arreglarse desde la lectura de la base, ver la opción na.strings de la función fread). Transformando esta variable a números, R genera NA (not available) en culaquier celda que tenga texto. Esta es la manera en la que R maneja los valores perdidos. Finalmente, se obtiene el número de valores perdidos y guardamos una base de datos nueva sin valores perdidos en la variable Magnitud.

DATOS <- DATOS0
DATOS$Magnitud <- as.numeric(DATOS0$Magnitud)
## Warning: NAs introduced by coercion
sum(is.na(DATOS$Magnitud))
## [1] 17649
DATOS <- DATOS[!is.na(DATOS$Magnitud),]

Manejo de datos y análisis exploratorio

Analizaremos las magnitudes por años, meses, etc. Por lo que es necesario crear algunas variables nuevas a partir de la fecha en la que se registró cada evento sísmico en la base de datos. Es importante mencionar que la función fred detectó correctamente que la variabe Fecha es de tipo date, que es el tipo de variable en el que R maneja a las fechas. Considerando lo anterior, se extrae la información de la siguiente manera.

# SE EXTRAE EL AÑO, MES (en español) Y SE GUARDAN EN VARIABLES NUEVAS 
DATOS$year <- year(DATOS$Fecha)
meses_esp <- c("Ene","Feb","Mar","Abr","May","Jun","Jul", "Ago","Sep","Oct","Nov","Dic")
DATOS$mes <- factor(meses_esp[month(DATOS$Fecha)], levels = meses_esp)

# SE CREA LA VARIABLE DIA_MES, SIN CONSIDERAR EL AÑO
lv <- NULL
for (j in 1:12) {
  lv <- c(lv, paste0(meses_esp[j], "_", formatC(1:31, width = 2, format = "d", flag = "0"))) 
}
DATOS$DIA_MES <- factor(paste0(DATOS$mes, "_", format(DATOS$Fecha, format = "%d")), levels = lv)

Se tiene información para 108 años de un total de 124, por lo que hay años en los que no se registró información. A continuación se grafica la magnitud de los sismos por año: mínima, promedio y máxima.

tb1 <- tapply(DATOS$Magnitud, DATOS$year, min)
tb2 <- tapply(DATOS$Magnitud, DATOS$year, mean)
tb3 <- tapply(DATOS$Magnitud, DATOS$year, max)

# GRÁFICA 
my <- as.numeric(names(tb1))
plot(my, tb1, type = "l", lwd = 2, axes = FALSE, 
     xlab = "Año", ylab = "Magnitud de los sismos", col = "blue", xlim = c(1900, 2025), ylim = c(0, 9))
lines(my, tb2, lwd = 2)
lines(my, tb3, lwd = 2, col = "red")
legend("bottomleft", c("min", "prom", "max"), 
       lwd = c(2, 2, 2), col = c("blue", "black", "red"), bty = "n")
axis(1, seq(1900, 2025, by = 5), las = 2, cex.axis = 0.8)
axis(2, seq(0, 10, length.out = 11), las = 2, cex.axis = 0.8)
box(lwd = 2)

Es posible observar directamente lo que se mencionaba anteriormente. La capacidad para medir los sismos ha ido evlolucionando gradualmente. Además, se observa que la década de 1970 fue un parteaguas. Hasta antes de esta década sólo se registraban sismos mayores a seis grados de magnitud. A mediados de los 70’s la capacidad aumentó y fue posible medir sismos de menores magnitudes, así como sismos fuertes con mayor precisión (se observa mayor varianza en los máximos).

Finalmente, daremos un vistazo a los sismos de magnitud máxima. Primero, porque es algo que puede ser de interés general y segundo, es la información que podría considerarse comparable desde 1900 hasta el periodo actual. Con este objetivo en mente, se construye una base de datos con toda la información de los sismos de magnitud máxima por año. Observando que si en una año determinado hay empates para el sismo de magnitud máxima, se incluye la información de los empates.

my_year <- unique(DATOS$year)
DATOSMAX <- NULL
for (i in my_year) {
  DATOSTMP <- DATOS[DATOS$year == i,]
  DATOSMAX <- rbind(DATOSMAX, DATOSTMP[which(DATOSTMP$Magnitud == max(DATOSTMP$Magnitud)), ])
}

Se realiza el histograma de las magnitudes máximas y se ajusta una normal con media 6.99 y desviación estándar 0.503.

x <- seq(5.5, 8.3, length.out = 500)
y <- dnorm(x, mean(DATOSMAX$Magnitud),
              sd(DATOSMAX$Magnitud))
hist(DATOSMAX$Magnitud, breaks = seq(5.5, 8.3, length.out = 12), 
     ylim = c(0, 0.8), main = "", probability = TRUE, border = "white", axes = FALSE, 
     xlab = "Magnitud máxima anual de los sismos en México desde 1900 a 2023", ylab = "Densidad estimada")
lines(x, y, col = "red", lwd = 2)
axis(1, seq(5.5, 8.3, length.out = 11))
axis(2, seq(0, 0.8, length.out = 11), las = 2)
box(lwd = 2)

La normal ajustada parece describir adecuadamente los datos. Sin embargo, los máximos por año muestran una liguera asimetría y la cola derecha es un poco más pesada que la izquierda.

Finalmente, se consideran los meses en los que ocurrieron los sismos máximos y se calculan los porcentajes correspondientes (se ordenan de mayor a menor).

tb2 <- table(DATOSMAX$mes)
pp <- sort(round(100*tb2/sum(tb2), 1), decreasing = TRUE)
rr <- barplot(pp, border = "white", las = 2, ylim = c(0, 16))
text(rr, pp+0.5, pp, cex = 0.8)
box(lwd = 2)

Diciembre, junio y septiembre son los meses en los que han ocurrido los sismos de magnitud más fuerte de cada año. Lo anterior considerando el periodo desde el año de 1900 a lo que va del 2023. El 100% son 116 sismos.