Se presenta el código fuente en R utilizado para la descarga, procesamiento y análisis de los datos sísmicos correspondientes a la Asesoría I:
############### Script Asesoria I ###############
# Librerias ----
library(readr)
library(dplyr)
library(stringr)
library(corrplot)
library(ggplot2)
library(maps)
library(leaflet)
library(htmltools)
library(ggcorrplot)
library(patchwork)
library(vcd)
library(FactoMineR)
library(factoextra)
library(ggrepel)
library(leaflet.extras)
library(sf)
library(crosstalk)
library(viridisLite)
library(tidyr)
# Obtención de datos ----
descargar_usgs_anual <- function(minlat, maxlat, minlon, maxlon,
anio_inicio, anio_fin,
minmag ) {
base_url <- "https://earthquake.usgs.gov/fdsnws/event/1/query"
lista_datos <- list()
for (anio in anio_inicio:anio_fin) {
starttime <- paste0(anio, "-01-01")
endtime <- paste0(anio + 1, "-01-01")
url <- paste0(
base_url,
"?format=csv",
"&starttime=", starttime,
"&endtime=", endtime,
"&minlatitude=", minlat,
"&maxlatitude=", maxlat,
"&minlongitude=", minlon,
"&maxlongitude=", maxlon,
"&minmagnitude=", minmag,
"&eventtype=earthquake",
"&orderby=time-asc"
)
cat("Descargando año:", anio, "\n")
datos_anio <- tryCatch(
read_csv(url, show_col_types = FALSE),
error = function(e) {
message("Error en año ", anio, ": ", e$message)
return(NULL)
}
)
if (!is.null(datos_anio) && nrow(datos_anio) > 0) {
lista_datos[[as.character(anio)]] <- datos_anio
}
}
bind_rows(lista_datos)
}
Base_CA <- descargar_usgs_anual(
minlat = 32,
maxlat = 42,
minlon = -125,
maxlon = -114,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 5
)
Base_CHILE <- descargar_usgs_anual(
minlat = -46,
maxlat = -17,
minlon = -76,
maxlon = -66,
anio_inicio = 2000,
anio_fin = 2025,
minmag = 5
)
# GUARDADO DE LA BASE CRUDA
getwd()
write.csv(Base_CHILE,"chile_limpia.csv", row.names = FALSE)
write.csv(Base_CA, "california_limpia.csv", row.names = FALSE)
# verificación de datos faltantes
variables_clave <- c(
"id", "time", "latitude", "longitude", "place",
"depth", "mag", "magType", "type"
)
eliminadas_california <- sum(!complete.cases(california[, variables_clave]))
eliminadas_chile <- sum(!complete.cases(chile[, variables_clave]))
eliminadas_california #Sin obs con datos faltantes
eliminadas_chile #Sin obs con datos faltantes
# duplicados
duplicados_exactos_california <- california[duplicated(california), ]
nrow(duplicados_exactos_california)
duplicados_exactos_chile <- chile[duplicated(chile), ]
nrow(duplicados_exactos_chile)
# Transformacion Escalas magnitud ----
table(chile$magType)
table(california$magType)
california %>%
count(magType) %>%
mutate(porcentaje = round(100 * n / sum(n), 2))
chile %>%
count(magType) %>%
mutate(porcentaje = round(100 * n / sum(n), 2))
## Homogenizacion Chile ----
chile <- chile %>%
mutate(
Mw_hom = case_when(
magType %in% c("mw", "mwb", "mwc", "mwr", "mww") ~ mag,
magType == "ml" & depth <= 50 ~ 0.80 * mag + 1.15,
magType == "ml" & depth > 50 ~ 0.94 * mag + 0.30,
magType == "ms" ~ 0.74 * mag + 1.60,
magType == "mb" ~ 1.04 * mag - 0.02,
magType %in% c("m", "md", "mc") ~ mag,
TRUE ~ NA_real_
)
)
## Homogenizacion California -- NO REALIZADA -- Filtracion de datos solo a Mw
california <- california[
tolower(california$magType) %in% c("mw","mwr"),
]
table(california$magType)
# Categorización Magnitud y Profundidad ----
## Magnitud ----
chile <- chile %>% mutate(mag_class=ifelse(Mw_hom < 4.0,"Menor",
ifelse(Mw_hom < 5.0, "Ligera",
ifelse(Mw_hom < 6.0,"Moderada",
ifelse(Mw_hom < 7.0, "Fuerte",
ifelse(Mw_hom < 8.0, "Mayor","Gran Terremoto"))))))
california <- california %>% mutate(mag_class=ifelse(mag < 4.0,"Menor",
ifelse(mag < 5.0, "Ligera",
ifelse(mag < 6.0,"Moderada",
ifelse(mag < 7.0, "Fuerte",
ifelse(mag < 8.0, "Mayor","Gran Terremoto"))))))
table(chile$mag_class)
table(california$mag_class)
## Profundidad ----
chile <- chile %>% mutate(depth_class=ifelse(depth <= 20, "Superficial Superior",
ifelse(depth <= 70, "Superficial Inferior",
ifelse(depth <= 300, "Intermedio", "Profundo"))))
california <- california %>% mutate(depth_class=ifelse(depth <= 20, "Superficial Superior",
ifelse(depth <= 70, "Superficial Inferior",
ifelse(depth <= 300, "Intermedio", "Profundo"))))
table(chile$depth_class)
table(california$depth_class)
# Analisis Univariado ----
## Comparación Magnitudes ----
### Tablas ----
summary(california$mag)
summary(chile$Mw_hom)
sd(california$mag)
sd(chile$Mw_hom)
quantile(california$mag,probs = c(0.25, 0.50, 0.75, 0.90,0.95), na.rm = TRUE,type = 7)
quantile(chile$Mw_hom,probs = c(0.25, 0.50, 0.75, 0.90,0.95), na.rm = TRUE,type = 7)
### Histogramas ----
ancho_bin <- 0.25
plot_ca <- ggplot(california, aes(x = mag)) +
geom_histogram(binwidth = ancho_bin, color = "white", fill = "#E69F00", alpha = 0.7) +
geom_density(aes(y = after_stat(count * ancho_bin)),
color = "#b37b00", linewidth = 1.2, fill = "#E69F00", alpha = 0.2) +
labs(
title = "California",
x = "Magnitud",
y = "Número de Sismos"
) +
theme_minimal(base_size = 12)
plot_chi <- ggplot(chile, aes(x = Mw_hom)) +
geom_histogram(binwidth = ancho_bin, color = "white", fill = "#0072B2", alpha = 0.7) +
geom_density(aes(y = after_stat(count * ancho_bin)),
color = "#00558a", linewidth = 1.2, fill = "#0072B2", alpha = 0.2) +
labs(
title = "Chile",
x = "Magnitud",
y = "Número de Sismos"
) +
theme_minimal(base_size = 12)
plot_final <- plot_ca + plot_chi +
plot_annotation(
title = "Distribución de Magnitud",
subtitle = "Número de sismos según Magnitud entre California y Chile",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 12, hjust = 0.5)
)
)
plot_final
### BoxPlots ----
df_california <- california %>%
select(Magnitud = mag) %>%
mutate(Region = "California")
df_chile <- chile %>%
select(Magnitud = Mw_hom) %>%
mutate(Region = "Chile")
df_terremotos <- bind_rows(df_california, df_chile)
ggplot(df_terremotos, aes(x = Region, y = Magnitud, fill = Region)) +
geom_boxplot(
alpha = 0.7,
width = 0.5,
color = "black",
outlier.color = "red",
outlier.alpha = 0.5,
outlier.size = 2
) +
scale_fill_manual(values = c("California" = "#E69F00", "Chile" = "#0072B2")) +
labs(
title = "Comparación de Magnitudes Sísmicas",
subtitle = "Distribución de Magnitudes mediante BoxPlots",
x = NULL,
y = "Magnitud"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(size = 14, face = "bold", color = "black")
)
## Comparación Profundidad ----
### Tablas ----
summary(california$depth)
summary(chile$depth)
sd(california$depth)
sd(chile$depth)
quantile(california$depth,probs = c(0.25, 0.50, 0.75, 0.90,0.95), na.rm = TRUE,type = 7)
quantile(chile$depth,probs = c(0.25, 0.50, 0.75, 0.90,0.95), na.rm = TRUE,type = 7)
### Histogramas ----
ancho_bin <- 4.5
plot_ca <- ggplot(california, aes(x = depth)) +
geom_histogram(binwidth = ancho_bin, color = "white", fill = "#E69F00", alpha = 0.7) +
geom_density(aes(y = after_stat(count * ancho_bin)),
color = "#b37b00", linewidth = 1.2, fill = "#E69F00", alpha = 0.2) +
labs(
title = "California",
x = "Profundidad",
y = "Número de Sismos"
) +
theme_minimal(base_size = 12)
plot_chi <- ggplot(chile, aes(x = depth)) +
geom_histogram(binwidth = ancho_bin, color = "white", fill = "#0072B2", alpha = 0.7) +
geom_density(aes(y = after_stat(count * ancho_bin)),
color = "#00558a", linewidth = 1.2, fill = "#0072B2", alpha = 0.2) +
labs(
title = "Chile",
x = "Profundidad",
y = "Número de Sismos"
) +
theme_minimal(base_size = 12)
plot_final <- plot_ca + plot_chi +
plot_annotation(
title = "Distribución de Profundidad",
subtitle = "Número de sismos según Profundidad entre California y Chile",
theme = theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 12, hjust = 0.5)
)
)
plot_final
### BoxPlots ----
df_california <- california %>%
select(Profundidad = depth) %>%
mutate(Region = "California")
df_chile <- chile %>%
select(Profundidad = depth) %>%
mutate(Region = "Chile")
df_terremotos <- bind_rows(df_california, df_chile)
ggplot(df_terremotos, aes(x = Region, y = Profundidad, fill = Region)) +
geom_boxplot(
alpha = 0.7,
width = 0.5,
color = "black",
outlier.color = "red",
outlier.alpha = 0.5,
outlier.size = 2
) +
scale_fill_manual(values = c("California" = "#E69F00", "Chile" = "#0072B2")) +
labs(
title = "Comparación de Profundidades Sísmicas",
subtitle = "Distribución de cuartiles, mediana y valores atípicos",
x = NULL,
y = "Profundidad (km)"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(size = 14, face = "bold", color = "black")
)
# Analisis Bivariado ----
## California ----
### Asociación Magnitud-Profundidad-Latitud-Longitud ----
plot_corr <- function(base, metodo, titulo) {
var <- base[, c("mag", "depth", "latitude", "longitude")]
colnames(var) <- c("Magnitud", "Profundidad", "Latitud", "Longitud")
mat <- cor(var, use = "complete.obs", method = metodo)
ggcorrplot(
mat,
method = "square",
type = "lower",
lab = TRUE,
lab_size = 4,
digits = 2,
colors = c("#B2182B", "white", "#2166AC"),
outline.color = "white",
ggtheme = theme_minimal(base_size = 11),
title = titulo,
legend.title = "r"
) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
legend.position = "right"
)
}
ca_pearson <- plot_corr(california, "pearson", "Pearson")
ca_spearman <- plot_corr(california, "spearman", "Spearman")
ca_kendall <- plot_corr(california, "kendall", "Kendall")
figura_california <- (ca_pearson | ca_spearman | ca_kendall) +
plot_annotation(
title = "Matrices de correlación — California",
subtitle = "Variables: Magnitud, Profundidad, Latitud, Longitud",
theme = theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 13),
plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray40")
)
)
figura_california
## Chile ----
### Asociación Magnitud-Profundidad-Latitud-Longitud ----
chile_aux <- chile
chile_aux$mag <- chile_aux$Mw_hom
ch_pearson <- plot_corr(chile_aux, "pearson", "Pearson")
ch_spearman <- plot_corr(chile_aux, "spearman", "Spearman")
ch_kendall <- plot_corr(chile_aux, "kendall", "Kendall")
figura_chile <- (ch_pearson | ch_spearman | ch_kendall) +
plot_annotation(
title = "Matrices de correlación — Chile",
subtitle = "Variables: Magnitud, Profundidad, Latitud, Longitud",
theme = theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 13),
plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray40")
)
)
figura_chile
# Analisis exploratorio ----
## Mapas Georreferenciados ----
### Mapas Interactivos ----
# Placas tectónicas
# Primero se intenta cargar GPlates. Si el servidor no responde,
# se usa un trazado local aproximado para que las placas SIGAN apareciendo.
crear_placas_locales <- function() {
margen_chile <- st_linestring(matrix(c(
-75.5, -46.0,
-74.8, -42.0,
-74.0, -38.0,
-73.2, -34.0,
-72.4, -30.0,
-71.4, -26.0,
-70.2, -22.0,
-69.0, -18.0
), ncol = 2, byrow = TRUE))
falla_san_andres <- st_linestring(matrix(c(
-124.0, 40.2,
-123.0, 39.0,
-122.3, 38.0,
-121.4, 37.0,
-120.5, 36.0,
-119.6, 35.2,
-118.2, 34.4,
-116.7, 33.3,
-115.4, 32.6
), ncol = 2, byrow = TRUE))
zona_mendocino <- st_linestring(matrix(c(
-127.0, 40.3,
-125.0, 40.25,
-123.0, 40.15
), ncol = 2, byrow = TRUE))
st_sf(
nombre = c(
"Margen de subducción de Chile",
"Falla de San Andrés / límite transformante",
"Zona de Mendocino"
),
tipo = "Límite de placa tectónica",
fuente = "Trazado local aproximado usado como respaldo cuando GPlates no responde",
geometry = st_sfc(margen_chile, falla_san_andres, zona_mendocino, crs = 4326)
)
}
placas_actuales <- tryCatch(
{
st_read(
"https://gws.gplates.org/topology/plate_boundaries?time=0&model=ZAHIROVIC2022",
quiet = TRUE
) %>%
st_transform(4326)
},
error = function(e) {
crear_placas_locales()
}
)
# Seguridad extra: si GPlates devuelve algo vacío, se usan las placas locales.
if (nrow(placas_actuales) == 0) {
placas_actuales <- crear_placas_locales()
}
# Leyenda continua para el mapa de calor
leyenda_heatmap <- function(titulo = "Mapa de calor") {
HTML(paste0(
"<div style='background:white; padding:8px; border-radius:5px;",
"box-shadow:0 0 5px rgba(0,0,0,0.3); font-size:12px;'>",
"<b>", titulo, "</b><br>",
"<div style='width:170px; height:14px; margin-top:6px; margin-bottom:4px;",
"background: linear-gradient(to right, blue, cyan, green, yellow, red);",
"border:1px solid #999;'></div>",
"<div style='display:flex; justify-content:space-between; width:170px;'>",
"<span>Baja</span><span>Alta</span>",
"</div>",
"<span>Concentración / intensidad sísmica</span>",
"</div>"
))
}
#### MAPA CHILE CON FILTROS + PLACAS + EPICENTROS + MAPA DE CALOR ----
# Convertir fecha
chile$time <- as.POSIXct(
chile$time,
format = "%Y-%m-%dT%H:%M:%OSZ",
tz = "UTC"
)
# Año
chile$year <- as.integer(format(chile$time, "%Y"))
# Crear magnitud común para el mapa
# En Chile usamos magnitud homogenizada Mw_hom
chile$mag_mapa <- chile$Mw_hom
# Filtrar datos válidos
chile_mapa <- chile %>%
filter(
!is.na(longitude),
!is.na(latitude),
!is.na(depth),
!is.na(year),
!is.na(mag_mapa)
)
# ID único para filtros
chile_mapa$id_filtro <- paste0("chile_", seq_len(nrow(chile_mapa)))
# Popup
chile_mapa$popup_chile <- paste0(
"<b>Lugar:</b> ", chile_mapa$place, "<br>",
"<b>Fecha:</b> ", chile_mapa$time, "<br>",
"<b>Año:</b> ", chile_mapa$year, "<br>",
"<b>Magnitud Mw:</b> ", round(chile_mapa$mag_mapa, 2), "<br>",
"<b>Profundidad:</b> ", round(chile_mapa$depth, 2), " km"
)
# Placas tectónicas Chile
placas_chile <- st_crop(
placas_actuales,
xmin = -80,
xmax = -65,
ymin = -56,
ymax = -15
)
# Paleta de color por profundidad
pal_depth_chile <- colorNumeric(
palette = viridis(256),
domain = chile_mapa$depth
)
# Nombres de placas para Chile
nombres_placas_chile <- data.frame(
nombre = c("Placa de Nazca", "Placa Sudamericana"),
lon = c(-78.5, -68.5),
lat = c(-31, -31)
)
# Leyenda manual
leyenda_mag_chile <- HTML("
<div style='background:white; padding:8px; border-radius:5px;
box-shadow:0 0 5px rgba(0,0,0,0.3); font-size:12px;'>
<b>Tamaño del punto</b><br>
<span style='display:inline-block; width:8px; height:8px;
border-radius:50%; background:gray;'></span> Menor magnitud<br>
<span style='display:inline-block; width:14px; height:14px;
border-radius:50%; background:gray;'></span> Magnitud media<br>
<span style='display:inline-block; width:22px; height:22px;
border-radius:50%; background:gray;'></span> Mayor magnitud
</div>
")
# SharedData para filtros
chile_sd <- SharedData$new(
chile_mapa,
key = ~id_filtro,
group = "sismos_chile"
)
# Mapa Chile
mapa_chile_leaflet <- bscols(
widths = c(2, 10),
list(
h4("Filtros Chile"),
filter_slider(
id = "filtro_anio_chile",
label = "Año",
sharedData = chile_sd,
column = ~year,
step = 1,
width = "100%"
),
filter_slider(
id = "filtro_mag_chile",
label = "Magnitud Mw",
sharedData = chile_sd,
column = ~mag_mapa,
step = 0.1,
width = "100%"
),
filter_slider(
id = "filtro_depth_chile",
label = "Profundidad (km)",
sharedData = chile_sd,
column = ~depth,
step = 1,
width = "100%"
)
),
leaflet(chile_sd, height = "780px") %>%
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa claro"
) %>%
addProviderTiles(
providers$CartoDB.DarkMatter,
group = "Mapa oscuro"
) %>%
addPolylines(
data = placas_chile,
color = "black",
weight = 6,
opacity = 0.9,
group = "Placas tectónicas"
) %>%
addPolylines(
data = placas_chile,
color = "#FF10F0",
weight = 3,
opacity = 1,
popup = "Límite de placa tectónica - GPlates ZAHIROVIC2022",
group = "Placas tectónicas"
) %>%
addLabelOnlyMarkers(
data = nombres_placas_chile,
lng = ~lon,
lat = ~lat,
label = ~nombre,
labelOptions = labelOptions(
noHide = TRUE,
direction = "center",
textOnly = TRUE,
style = list(
"font-weight" = "bold",
"font-size" = "15px",
"color" = "black",
"text-shadow" = "1px 1px 2px white"
)
),
group = "Nombres placas"
) %>%
addHeatmap(
data = chile_mapa,
lng = ~longitude,
lat = ~latitude,
intensity = ~mag_mapa,
radius = 15,
blur = 10,
max = max(chile_mapa$mag_mapa, na.rm = TRUE),
minOpacity = 0.65,
gradient = c(
"0.10" = "blue",
"0.30" = "cyan",
"0.50" = "green",
"0.70" = "yellow",
"1.00" = "red"
),
group = "Mapa de calor"
) %>%
addCircleMarkers(
lng = ~longitude,
lat = ~latitude,
radius = ~pmax((mag_mapa - 4) * 4, 3),
color = "black",
stroke = TRUE,
weight = 0.7,
opacity = 0.9,
fillColor = ~pal_depth_chile(depth),
fillOpacity = 0.75,
popup = ~popup_chile,
group = "Epicentros"
) %>%
addLegend(
position = "bottomright",
pal = pal_depth_chile,
values = chile_mapa$depth,
title = "Profundidad (km)",
opacity = 1
) %>%
addControl(
leyenda_mag_chile,
position = "bottomleft"
) %>%
addControl(
leyenda_heatmap("Mapa de calor"),
position = "bottomleft"
) %>%
addLayersControl(
position = "topright",
baseGroups = c("Mapa claro", "Mapa oscuro"),
overlayGroups = c(
"Placas tectónicas",
"Nombres placas",
"Epicentros",
"Mapa de calor"
),
options = layersControlOptions(collapsed = TRUE)
) %>%
hideGroup("Mapa de calor") %>%
fitBounds(
lng1 = -76,
lat1 = -46,
lng2 = -66,
lat2 = -17
)
)
mapa_chile_leaflet
#### MAPA CALIFORNIA CON FILTROS + PLACAS + EPICENTROS + MAPA DE CALOR ----
# Convertir fecha
california$time <- as.POSIXct(
california$time,
format = "%Y-%m-%dT%H:%M:%OSZ",
tz = "UTC"
)
# Año
california$year <- as.integer(format(california$time, "%Y"))
# Crear magnitud común para el mapa
# En California usamos mag, porque ya fue filtrada a Mw
california$mag_mapa <- california$mag
# Filtrar datos válidos
california_mapa <- california %>%
filter(
!is.na(longitude),
!is.na(latitude),
!is.na(depth),
!is.na(year),
!is.na(mag_mapa)
)
# ID único para filtros
california_mapa$id_filtro <- paste0("ca_", seq_len(nrow(california_mapa)))
# Popup
california_mapa$popup_ca <- paste0(
"<b>Lugar:</b> ", california_mapa$place, "<br>",
"<b>Fecha:</b> ", california_mapa$time, "<br>",
"<b>Año:</b> ", california_mapa$year, "<br>",
"<b>Magnitud Mw:</b> ", round(california_mapa$mag_mapa, 2), "<br>",
"<b>Profundidad:</b> ", round(california_mapa$depth, 2), " km"
)
# Placas tectónicas California
placas_california <- st_crop(
placas_actuales,
xmin = -127,
xmax = -110,
ymin = 30,
ymax = 43
)
# Paleta de color por profundidad
pal_depth_ca <- colorNumeric(
palette = viridis(256),
domain = california_mapa$depth
)
# Nombres de placas para California
nombres_placas_ca <- data.frame(
nombre = c("Placa del Pacífico", "Placa Norteamericana"),
lon = c(-124.2, -115.5),
lat = c(36.5, 36.5)
)
# Leyenda manual California
leyenda_mag_ca <- HTML("
<div style='background:white; padding:8px; border-radius:5px;
box-shadow:0 0 5px rgba(0,0,0,0.3); font-size:12px;'>
<b>Tamaño del punto</b><br>
<span style='display:inline-block; width:8px; height:8px;
border-radius:50%; background:gray;'></span> Menor magnitud<br>
<span style='display:inline-block; width:14px; height:14px;
border-radius:50%; background:gray;'></span> Magnitud media<br>
<span style='display:inline-block; width:22px; height:22px;
border-radius:50%; background:gray;'></span> Mayor magnitud
</div>
")
# SharedData
california_sd <- SharedData$new(
california_mapa,
key = ~id_filtro,
group = "sismos_california"
)
# Mapa California
mapa_california_leaflet <- bscols(
widths = c(2, 10),
list(
h4("Filtros California"),
filter_slider(
id = "filtro_anio_ca",
label = "Año",
sharedData = california_sd,
column = ~year,
step = 1,
width = "100%"
),
filter_slider(
id = "filtro_mag_ca",
label = "Magnitud Mw",
sharedData = california_sd,
column = ~mag_mapa,
step = 0.1,
width = "100%"
),
filter_slider(
id = "filtro_depth_ca",
label = "Profundidad (km)",
sharedData = california_sd,
column = ~depth,
step = 1,
width = "100%"
)
),
leaflet(california_sd, height = "650px") %>%
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa claro"
) %>%
addProviderTiles(
providers$CartoDB.DarkMatter,
group = "Mapa oscuro"
) %>%
addPolylines(
data = placas_california,
color = "black",
weight = 6,
opacity = 0.9,
group = "Placas tectónicas"
) %>%
addPolylines(
data = placas_california,
color = "#FF10F0",
weight = 3,
opacity = 1,
popup = "Límite de placa tectónica - GPlates ZAHIROVIC2022",
group = "Placas tectónicas"
) %>%
addLabelOnlyMarkers(
data = nombres_placas_ca,
lng = ~lon,
lat = ~lat,
label = ~nombre,
labelOptions = labelOptions(
noHide = TRUE,
direction = "center",
textOnly = TRUE,
style = list(
"font-weight" = "bold",
"font-size" = "15px",
"color" = "black",
"text-shadow" = "1px 1px 2px white"
)
),
group = "Nombres placas"
) %>%
addHeatmap(
data = california_mapa,
lng = ~longitude,
lat = ~latitude,
intensity = ~mag_mapa,
radius = 22,
blur = 12,
max = max(california_mapa$mag_mapa, na.rm = TRUE),
minOpacity = 0.70,
gradient = c(
"0.10" = "blue",
"0.30" = "cyan",
"0.50" = "green",
"0.70" = "yellow",
"1.00" = "red"
),
group = "Mapa de calor"
) %>%
addCircleMarkers(
lng = ~longitude,
lat = ~latitude,
radius = ~pmax((mag_mapa - 4) * 4, 3),
color = "black",
stroke = TRUE,
weight = 0.7,
opacity = 0.9,
fillColor = ~pal_depth_ca(depth),
fillOpacity = 0.75,
popup = ~popup_ca,
group = "Epicentros"
) %>%
addLegend(
position = "bottomright",
pal = pal_depth_ca,
values = california_mapa$depth,
title = "Profundidad (km)",
opacity = 1
) %>%
addControl(
leyenda_mag_ca,
position = "bottomleft"
) %>%
addControl(
leyenda_heatmap("Mapa de calor"),
position = "bottomleft"
) %>%
addLayersControl(
position = "topright",
baseGroups = c("Mapa claro", "Mapa oscuro"),
overlayGroups = c(
"Placas tectónicas",
"Nombres placas",
"Epicentros",
"Mapa de calor"
),
options = layersControlOptions(collapsed = TRUE)
) %>%
hideGroup("Mapa de calor") %>%
fitBounds(
lng1 = -125,
lat1 = 32,
lng2 = -114,
lat2 = 42
)
)
mapa_california_leaflet
### Analisis Sismos Intermedios: Chile ----
sismos_rango_chile <- chile %>%
filter(depth >= 95 & depth <= 125)
nrow(sismos_rango_chile)
#Mapeo Georreferenciado
mapa_chile_rango <- map_data("world") %>% filter(region == "Chile")
ggplot() +
geom_polygon(data = mapa_chile_rango, aes(x = long, y = lat, group = group),
fill = "gray93", color = "gray50", linewidth = 0.4) +
geom_point(data = sismos_rango_chile, aes(x = longitude, y = latitude),
color = "darkblue", size = 1, alpha = 0.5) +
coord_fixed(ratio = 1.3, xlim = c(-76, -66), ylim = c(-56, -17)) +
labs(
title = "Sismos en el Rango de 95-125 km de profundidad: Chile",
subtitle = "Dispersión geográfica de sismicidad de profundidad intermedia",
x = "Longitud",
y = "Latitud"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5),
panel.grid.major = element_line(color = "gray90")
)
### Analisis Sismos Superficial Superior: California ----
sismos_rango_california <- california %>%
filter(depth >= 25 & depth <= 30)
nrow(sismos_rango_california)
# Mapeo Georreferenciado
mapa_california_rango <- map_data("state") %>% filter(region == "california")
ggplot() +
geom_polygon(data = mapa_california_rango, aes(x = long, y = lat, group = group),
fill = "gray93", color = "gray50", linewidth = 0.4) +
geom_point(data = sismos_rango_california, aes(x = longitude, y = latitude),
color = "darkblue", size = 1, alpha = 0.5) +
coord_fixed(ratio = 1.3, xlim = c(-125, -114), ylim = c(32, 42)) +
labs(
title = "Sismos en el Rango de 25-30 km de profundidad: California",
subtitle = "Dispersión Geográfica de sismicidad de profundidad superficial superior",
x = "Longitud",
y = "Latitud"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 11, hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5),
panel.grid.major = element_line(color = "gray90")
)
## Series de Tiempo ----
### Construccion Serie ----
construir_serie_mensual <- function(df, categoria, anio_inicio = 2000, anio_fin = 2025) {
sub <- df %>% filter(mag_class == categoria)
sub$year_month <- format(sub$time, "%Y-%m")
todos_los_meses <- format(
seq(as.Date(paste0(anio_inicio, "-01-01")),
as.Date(paste0(anio_fin, "-12-01")),
by = "month"),
"%Y-%m"
)
conteo <- sub %>%
count(year_month) %>%
right_join(data.frame(year_month = todos_los_meses), by = "year_month") %>%
mutate(n = ifelse(is.na(n), 0, n)) %>%
arrange(year_month)
conteo$categoria <- categoria
conteo
}
categorias <- c("Menor", "Ligera", "Moderada", "Fuerte", "Mayor", "Gran Terremoto")
series_california <- bind_rows(lapply(categorias, function(cat) construir_serie_mensual(california, cat)))
series_california$pais <- "California"
series_chile <- bind_rows(lapply(categorias, function(cat) construir_serie_mensual(chile, cat)))
series_chile$pais <- "Chile"
series_todas <- bind_rows(series_california, series_chile)
# Fijar el orden correcto de las categorías
series_todas <- series_todas %>%
mutate(
categoria = case_when(
tolower(categoria) == "menor" ~ "Menor",
tolower(categoria) == "ligera" ~ "Ligera",
tolower(categoria) == "moderada" ~ "Moderada",
tolower(categoria) == "fuerte" ~ "Fuerte",
tolower(categoria) == "mayor" ~ "Mayor",
tolower(categoria) == "gran terremoto" ~ "Gran Terremoto",
TRUE ~ as.character(categoria)
),
categoria = factor(categoria, levels = categorias)
)
series_todas$year_month_date <- as.Date(paste0(series_todas$year_month, "-01"))
sum(is.na(series_todas$categoria))
### Funcion Breaks adaptos para eje y ----
breaks_enteros_adaptados <- function(x) {
rango <- ceiling(x[2]) - floor(x[1])
if (rango <= 1) {
seq(0, 3, by = 1)
} else if (rango <= 10) {
seq(floor(x[1]), ceiling(x[2]), by = 1)
} else if (rango <= 30) {
seq(floor(x[1]), ceiling(x[2]), by = 5)
} else if (rango <= 60) {
seq(floor(x[1]), ceiling(x[2]), by = 10)
} else {
seq(floor(x[1]), ceiling(x[2]), by = 20)
}
}
### Gráfico comparativo final: California vs. Chile ----
categorias_incluidas <- c("Moderada", "Fuerte", "Mayor", "Gran Terremoto")
rangos_magnitud <- c(
"Menor" = "(<4.0)",
"Ligera" = "(4.0-4.9)",
"Moderada" = "(5.0-5.9)",
"Fuerte" = "(6.0-6.9)",
"Mayor" = "(7.0-7.9)",
"Gran Terremoto" = "(\u22658.0)"
)
series_todas$categoria_con_rango <- paste(series_todas$categoria,
rangos_magnitud[as.character(series_todas$categoria)])
series_todas$etiqueta_panel <- paste(series_todas$categoria_con_rango, "-", series_todas$pais)
categorias_con_rango_orden <- paste(categorias_incluidas, rangos_magnitud[categorias_incluidas])
orden_paneles <- as.vector(t(outer(categorias_con_rango_orden, c("California", "Chile"),
FUN = function(c, p) paste(c, "-", p))))
series_todas$etiqueta_panel <- factor(series_todas$etiqueta_panel, levels = orden_paneles)
ggplot(series_todas %>% filter(categoria %in% categorias_incluidas),
aes(x = year_month_date, y = n, color = pais)) +
geom_line() +
facet_wrap(~ etiqueta_panel, scales = "free_y", ncol = 2) +
scale_color_manual(values = c("California" = "steelblue", "Chile" = "firebrick")) +
scale_x_date(date_breaks = "6 months", date_labels = "%b %Y",
limits = as.Date(c("2000-01-01", "2025-12-01"))) +
scale_y_continuous(
breaks = breaks_enteros_adaptados,
limits = function(x) c(0, ifelse(ceiling(x[2]) <= 1, 3, ceiling(x[2]))),
expand = expansion(mult = c(0, 0.05))
) +
labs(
title = "Sismos por mes, California vs. Chile, según categoría de magnitud",
x = "Fecha", y = "N° de sismos en el mes"
) +
theme_minimal() +
theme(
strip.text = element_text(size = 7),
axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1, size = 5),
plot.margin = margin(t = 10, r = 10, b = 20, l = 10),
legend.position = "none"
)
### Descomposicion Clasica Aditiva ----
#### Chile - Moderada ----
serie_chile_moderada <- series_todas %>%
filter(pais == "Chile", categoria == "Moderada") %>%
arrange(year_month_date)
ts_chile_moderada <- ts(serie_chile_moderada$n, start = c(2000, 1), frequency = 12)
desc_chile_moderada <- decompose(ts_chile_moderada, type = "additive")
plot(desc_chile_moderada)
mtext("Chile - Sismos moderados (5.0 - 5.9)", side = 3, line = -3, outer = TRUE, font = 2)
#### Chile - Fuerte ----
serie_chile_fuerte <- series_todas %>%
filter(pais == "Chile", categoria == "Fuerte") %>%
arrange(year_month_date)
ts_chile_fuerte <- ts(serie_chile_fuerte$n, start = c(2000, 1), frequency = 12)
desc_chile_fuerte <- decompose(ts_chile_fuerte, type = "additive")
plot(desc_chile_fuerte)
mtext("Chile - Sismos fuertes (6.0 - 6.9)", side = 3, line = -3, outer = TRUE, font = 2)
#### California - Moderada ----
serie_california_moderada <- series_todas %>%
filter(pais == "California", categoria == "Moderada") %>%
arrange(year_month_date)
ts_california_moderada <- ts(serie_california_moderada$n, start = c(2000, 1), frequency = 12)
desc_california_moderada <- decompose(ts_california_moderada, type = "additive")
plot(desc_california_moderada)
mtext("California - Sismos moderados (5.0 - 5.9)", side = 3, line = -3, outer = TRUE, font = 2)
## Conteos Normalizados ----
### Parametros Generales ----
zonas <- c("Chile", "California")
niveles_magnitud <- c(
"5.0 - 5.9",
"6.0 - 6.9",
"7.0 - 7.9",
"8.0 - 8.9",
"9.0 - 9.9",
"10.0 o más"
)
niveles_profundidad <- c(
"0 - 20 km",
"20 - 70 km",
"70 - 300 km",
"300 km o más"
)
# Coordenadas de captura usadas en el trabajo.
# Chile: Oeste -76, Este -66, Sur -46, Norte -17
# California: Oeste -125, Este -114, Sur 32, Norte 42
crear_area_bbox <- function(oeste, este, sur, norte) {
poligono <- st_polygon(list(matrix(c(
oeste, sur,
este, sur,
este, norte,
oeste, norte,
oeste, sur
), ncol = 2, byrow = TRUE)))
sf_poligono <- st_sfc(poligono, crs = 4326)
as.numeric(st_area(st_transform(sf_poligono, 6933))) / 10^6
}
areas_zonas <- tibble(
zona = zonas,
area_km2 = c(
crear_area_bbox(oeste = -76, este = -66, sur = -46, norte = -17),
crear_area_bbox(oeste = -125, este = -114, sur = 32, norte = 42)
)
)
### Funciones Auxiliares ----
clasificar_magnitud <- function(magnitud) {
case_when(
magnitud >= 5.0 & magnitud < 6.0 ~ "5.0 - 5.9",
magnitud >= 6.0 & magnitud < 7.0 ~ "6.0 - 6.9",
magnitud >= 7.0 & magnitud < 8.0 ~ "7.0 - 7.9",
magnitud >= 8.0 & magnitud < 9.0 ~ "8.0 - 8.9",
magnitud >= 9.0 & magnitud < 10.0 ~ "9.0 - 9.9",
magnitud >= 10.0 ~ "10.0 o más",
TRUE ~ NA_character_
)
}
clasificar_profundidad <- function(profundidad) {
case_when(
profundidad >= 0 & profundidad < 20 ~ "0 - 20 km",
profundidad >= 20 & profundidad < 70 ~ "20 - 70 km",
profundidad >= 70 & profundidad < 300 ~ "70 - 300 km",
profundidad >= 300 ~ "300 km o más",
TRUE ~ NA_character_
)
}
preparar_base_zona <- function(base, nombre_zona, variable_mag) {
base %>%
mutate(
zona = nombre_zona,
time = as.POSIXct(time, format = "%Y-%m-%dT%H:%M:%OSZ", tz = "UTC"),
magnitud_usada = {{ variable_mag }},
rango_magnitud = factor(
clasificar_magnitud(magnitud_usada),
levels = niveles_magnitud
),
rango_profundidad = factor(
clasificar_profundidad(depth),
levels = niveles_profundidad
)
) %>%
select(zona, time, magnitud_usada, depth, rango_magnitud, rango_profundidad)
}
calcular_conteos_normalizados <- function(base, variables_grupo, combinaciones) {
base %>%
group_by(across(all_of(variables_grupo))) %>%
summarise(n_sismos = n(), .groups = "drop") %>%
right_join(combinaciones, by = variables_grupo) %>%
mutate(n_sismos = replace_na(n_sismos, 0)) %>%
group_by(zona, across(any_of("corte_magnitud"))) %>%
mutate(
total_zona_corte = sum(n_sismos),
proporcion = ifelse(total_zona_corte > 0, n_sismos / total_zona_corte, 0),
porcentaje = proporcion * 100
) %>%
ungroup() %>%
left_join(areas_zonas, by = "zona") %>%
mutate(
sismos_por_km2 = n_sismos / area_km2,
sismos_por_100mil_km2 = sismos_por_km2 * 100000
)
}
### Base comparativa común ----
base_comparativa <- bind_rows(
preparar_base_zona(chile, "Chile", Mw_hom),
preparar_base_zona(california, "California", mag)
)
### Conteos normalizados por categoría de magnitud ----
combinaciones_magnitud <- expand_grid(
zona = zonas,
rango_magnitud = factor(niveles_magnitud, levels = niveles_magnitud)
)
tabla_magnitud <- base_comparativa %>%
filter(!is.na(rango_magnitud)) %>%
calcular_conteos_normalizados(
variables_grupo = c("zona", "rango_magnitud"),
combinaciones = combinaciones_magnitud
) %>%
select(
zona,
rango_magnitud,
n_sismos,
total_zona = total_zona_corte,
proporcion,
porcentaje,
area_km2,
sismos_por_km2,
sismos_por_100mil_km2
) %>%
arrange(zona, rango_magnitud)
tabla_magnitud
# Gráfico 1: composición porcentual por magnitud
grafico_magnitud_porcentaje <- ggplot(
tabla_magnitud,
aes(x = rango_magnitud, y = porcentaje, fill = zona)
) +
geom_col(position = position_dodge(width = 0.9)) +
geom_text(
aes(label = paste0(round(porcentaje, 1), "%")),
position = position_dodge(width = 0.9),
vjust = -0.3,
size = 3.5
) +
scale_y_continuous(
limits = c(0, max(tabla_magnitud$porcentaje, na.rm = TRUE) * 1.15)
) +
labs(
title = "Composición porcentual de sismos por magnitud",
subtitle = "Control por tamaño muestral",
x = "Categoría de magnitud",
y = "Porcentaje dentro de cada zona (%)",
fill = "Zona"
) +
theme_minimal()
grafico_magnitud_porcentaje
# Gráfico 2: densidad espacial por magnitud
grafico_magnitud_area <- ggplot(
tabla_magnitud,
aes(x = rango_magnitud, y = sismos_por_100mil_km2, fill = zona)
) +
geom_col(position = position_dodge(width = 0.9)) +
geom_text(
aes(label = round(sismos_por_100mil_km2, 2)),
position = position_dodge(width = 0.9),
vjust = -0.3,
size = 3.5
) +
scale_y_continuous(
limits = c(0, max(tabla_magnitud$sismos_por_100mil_km2, na.rm = TRUE) * 1.15)
) +
labs(
title = "Densidad espacial de sismos por magnitud",
subtitle = "Sismos por cada 100.000 km²; mismo periodo de observación",
x = "Categoría de magnitud",
y = "Sismos por cada 100.000 km²",
fill = "Zona"
) +
theme_minimal()
grafico_magnitud_area
### Conteos normalizados por profundidad y corte de magnitud ----
cortes_magnitud <- c(5.0, 6.0)
base_profundidad_cortes <- bind_rows(
lapply(cortes_magnitud, function(corte) {
base_comparativa %>%
filter(!is.na(magnitud_usada), magnitud_usada >= corte) %>%
mutate(
corte_magnitud = factor(
paste0("Mw >= ", format(corte, nsmall = 1)),
levels = paste0("Mw >= ", format(cortes_magnitud, nsmall = 1))
)
)
})
) %>%
filter(!is.na(rango_profundidad))
combinaciones_profundidad_corte <- expand_grid(
zona = zonas,
corte_magnitud = factor(
paste0("Mw >= ", format(cortes_magnitud, nsmall = 1)),
levels = paste0("Mw >= ", format(cortes_magnitud, nsmall = 1))
),
rango_profundidad = factor(niveles_profundidad, levels = niveles_profundidad)
)
tabla_profundidad_cortes <- base_profundidad_cortes %>%
calcular_conteos_normalizados(
variables_grupo = c("zona", "corte_magnitud", "rango_profundidad"),
combinaciones = combinaciones_profundidad_corte
) %>%
select(
zona,
corte_magnitud,
rango_profundidad,
n_sismos,
total_zona_corte,
proporcion,
porcentaje,
area_km2,
sismos_por_km2,
sismos_por_100mil_km2
) %>%
arrange(zona, corte_magnitud, rango_profundidad)
tabla_profundidad_cortes
# Gráfico 3: composición porcentual por profundidad y corte de magnitud
grafico_profundidad_porcentaje <- ggplot(
tabla_profundidad_cortes,
aes(x = rango_profundidad, y = porcentaje, fill = corte_magnitud)
) +
geom_col(position = position_dodge(width = 0.9)) +
geom_text(
aes(label = paste0(round(porcentaje, 1), "%")),
position = position_dodge(width = 0.9),
vjust = -0.3,
size = 3.2
) +
facet_wrap(~ zona) +
scale_y_continuous(
limits = c(0, max(tabla_profundidad_cortes$porcentaje, na.rm = TRUE) * 1.15)
) +
labs(
title = "Composición porcentual de sismos por profundidad",
subtitle = "Control por tamaño muestral; comparación Mw >= 5.0 vs Mw >= 6.0",
x = "Categoría de profundidad",
y = "Porcentaje dentro de cada zona y corte (%)",
fill = "Corte de magnitud"
) +
theme_minimal()
grafico_profundidad_porcentaje
# Gráfico 4: densidad espacial por profundidad y corte de magnitud
grafico_profundidad_area <- ggplot(
tabla_profundidad_cortes,
aes(x = rango_profundidad, y = sismos_por_100mil_km2, fill = corte_magnitud)
) +
geom_col(position = position_dodge(width = 0.9)) +
geom_text(
aes(label = round(sismos_por_100mil_km2, 2)),
position = position_dodge(width = 0.9),
vjust = -0.3,
size = 3.2
) +
facet_wrap(~ zona) +
scale_y_continuous(
limits = c(0, max(tabla_profundidad_cortes$sismos_por_100mil_km2, na.rm = TRUE) * 1.15)
) +
labs(
title = "Densidad espacial de sismos por profundidad",
subtitle = "Sismos por cada 100.000 km²; comparación Mw >= 5.0 vs Mw >= 6.0",
x = "Categoría de profundidad",
y = "Sismos por cada 100.000 km²",
fill = "Corte de magnitud"
) +
theme_minimal()
grafico_profundidad_area
## Chi Cuadrado ----
# Agrupacion Magnitud Chile (3 niveles)
chile <- chile %>%
mutate(mag_class2 = factor(
case_when(
Mw_hom < 6.0 ~ "Moderada",
Mw_hom < 7.0 ~ "Fuerte",
TRUE ~ "Mayor o superior"
),
levels = c("Moderada", "Fuerte", "Mayor o superior")
))
# Magnitud California (2 niveles)
california <- california %>%
mutate(mag_class2 = factor(
case_when(
mag < 6.0 ~ "Moderada",
TRUE ~ "Fuerte o superior"
),
levels = c("Moderada", "Fuerte o superior")
))
# Profundidad a Factor Chile (3 niveles)
chile <- chile %>%
mutate(depth_factor = factor(
case_when(
depth <= 20 ~ "Superficial Superior",
depth <= 70 ~ "Superficial Inferior",
TRUE ~ "Intermedio"
),
levels = c("Superficial Superior", "Superficial Inferior", "Intermedio")
))
# Profundidad a Factor California (2 niveles)
california <- california %>%
mutate(depth_factor = factor(
case_when(
depth <= 20 ~ "Superficial Superior",
TRUE ~ "Superficial Inferior"
),
levels = c("Superficial Superior", "Superficial Inferior")
))
### Caso: Profundidad vs Magnitud - Chile ----
tabla1<-table("Profundidad" = chile$depth_factor, "Magnitud"= chile$mag_class2)
chisq.test(tabla1, simulate.p.value = TRUE)
assocstats(tabla1)$cramer
### Caso: Profundidad vs Magnitud - California ----
tabla2<-table("Profundidad" = california$depth_factor, "Magnitud"= california$mag_class2)
chisq.test(tabla2, simulate.p.value = TRUE)
assocstats(tabla2)$cramer
### Caso: Zona vs Magnitud ----
chile <- chile %>%
mutate(mag_class_conjunta = factor(
ifelse(Mw_hom < 6.0, "Moderada", "Fuerte o superior"),
levels = c("Moderada", "Fuerte o superior")
))
california <- california %>%
mutate(mag_class_conjunta = factor(
ifelse(mag < 6.0, "Moderada", "Fuerte o superior"),
levels = c("Moderada", "Fuerte o superior")
))
conjunto <- bind_rows(
chile %>% mutate(zona = "Chile") %>%
select(zona, mag_class_conjunta, depth_factor),
california %>% mutate(zona = "California") %>%
select(zona, mag_class_conjunta, depth_factor)
) %>%
mutate(zona = factor(zona, levels = c("Chile", "California")))
tabla3<-table("Zona" = conjunto$zona, "Magnitud" = conjunto$mag_class_conjunta)
chisq.test(tabla3, simulate.p.value = TRUE)
assocstats(tabla3)$cramer
### Caso: Zona VS Profundidad ----
tabla4<-table("Zona" = conjunto$zona, "Profundidad" = conjunto$depth_factor)
chisq.test(tabla4, simulate.p.value = TRUE)
assocstats(tabla4)$cramer
## Tablas de Contingencia ----
heatmap_tabla <- function(tabla, nombre_x, nombre_y, titulo, decimales = 0) {
dimnames(tabla) <- dimnames(tabla)
df <- as.data.frame(as.table(tabla))
colnames(df) <- c("Var1", "Var2", "Freq")
df <- df %>%
mutate(etiqueta = format(round(Freq, decimales), nsmall = decimales))
ggplot(df, aes(x = Var2, y = Var1, fill = Freq)) +
geom_tile(color = "white", linewidth = 0.8) +
geom_text(aes(label = etiqueta), size = 3.5, color = "black") +
scale_fill_gradient(low = "#EFF3FF", high = "#2171B5", name = "Frecuencia") +
labs(title = titulo, x = nombre_x, y = nombre_y) +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
axis.text.x = element_text(angle = 20, hjust = 1),
legend.position = "right"
)
}
esp_1 <- chisq.test(tabla1)$expected
dimnames(esp_1) <- dimnames(tabla1)
esp_2 <- chisq.test(tabla2)$expected
dimnames(esp_2) <- dimnames(tabla2)
esp_3 <- chisq.test(tabla3)$expected
dimnames(esp_3) <- dimnames(tabla3)
esp_4 <- chisq.test(tabla4)$expected
dimnames(esp_4) <- dimnames(tabla4)
p1_obs <- heatmap_tabla(tabla1, "Categoría de Magnitud", "Categoría de Profundidad",
"Caso 1: Observadas — Chile")
p1_esp <- heatmap_tabla(esp_1, "Categoría de Magnitud", "Categoría de Profundidad",
"Caso 1: Esperadas — Chile", decimales = 1)
fig_caso1 <- p1_obs | p1_esp
p2_obs <- heatmap_tabla(tabla2, "Categoría de Magnitud", "Categoría de Profundidad",
"Caso 2: Observadas — California")
p2_esp <- heatmap_tabla(esp_2, "Categoría de Magnitud", "Categoría de Profundidad",
"Caso 2: Esperadas — California", decimales = 1)
fig_caso2 <- p2_obs | p2_esp
p3_obs <- heatmap_tabla(tabla3, "Categoría de Magnitud", "Zona",
"Caso 3: Observadas")
p3_esp <- heatmap_tabla(esp_3, "Categoría de Magnitud", "Zona",
"Caso 3: Esperadas", decimales = 1)
fig_caso3 <- p3_obs | p3_esp
p4_obs <- heatmap_tabla(tabla4, "Categoría de Profundidad", "Zona",
"Caso 4: Observadas")
p4_esp <- heatmap_tabla(esp_4, "Categoría de Profundidad", "Zona",
"Caso 4: Esperadas", decimales = 1)
fig_caso4 <- p4_obs | p4_esp
fig_caso1
fig_caso2
fig_caso3
fig_caso4
## Analisis de Correspondencia ----
ac_caso4 <- CA(tabla4, graph = FALSE)
# Varianza explicada por dimensión
summary(ac_caso4)
ac_caso4$eig
coords_filas <- data.frame(
categoria = names(ac_caso4$row$coord),
Dim1 = as.numeric(ac_caso4$row$coord),
contrib = as.numeric(ac_caso4$row$contrib),
tipo = "Zona"
)
# Columnas (Profundidad) — matriz de 1 columna
coords_columnas <- data.frame(
categoria = rownames(ac_caso4$col$coord),
Dim1 = as.numeric(ac_caso4$col$coord),
contrib = as.numeric(ac_caso4$col$contrib),
tipo = "Profundidad"
)
coords_total <- rbind(coords_filas, coords_columnas)
# Verificar
print(coords_total)
var_dim1 <- round(ac_caso4$eig[1, 2], 1)
p_ac_clasico <- ggplot(coords_total,
aes(x = Dim1, y = 0,
color = tipo, shape = tipo)) +
geom_hline(yintercept = 0, linetype = "dashed",
color = "gray70", linewidth = 0.5) +
geom_vline(xintercept = 0, linetype = "dashed",
color = "gray70", linewidth = 0.5) +
geom_point(size = 5, alpha = 0.9) +
ggrepel::geom_text_repel(
aes(label = categoria),
size = 4,
fontface = "bold",
nudge_y = 0.05,
box.padding = 0.4,
show.legend = FALSE
) +
scale_color_manual(
values = c("Zona" = "#B2182B", "Profundidad" = "#2166AC"),
name = "Variable"
) +
scale_shape_manual(
values = c("Zona" = 17, "Profundidad" = 16),
name = "Variable"
) +
labs(
title = "Análisis de Correspondencia — Zona vs Profundidad",
subtitle = "Rojo (▲): Zona | Azul (●): Categoría de Profundidad",
x = paste0("Dimensión 1 (", var_dim1, "%)"),
y = "",
caption = ""
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 13),
plot.subtitle = element_text(hjust = 0.5, size = 10, color = "gray40"),
plot.caption = element_text(hjust = 0.5, size = 8, color = "gray50"),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
legend.position = "right"
) +
ylim(-0.2, 0.2)
p_ac_clasico
## PCA ----
# Chile
pca_chile <- prcomp(chile %>% select(Mw_hom, depth, latitude, longitude) %>% na.omit(),
center = TRUE, scale. = TRUE)
summary(pca_chile)
pca_chile$rotation
# California
pca_california <- prcomp(california %>% select(mag, depth, latitude, longitude) %>% na.omit(),
center = TRUE, scale. = TRUE)
summary(pca_california)
pca_california$rotation