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.
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:
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.
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.
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.
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.
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:
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.
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.
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),]
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.