ANEST
Cesar Jimenez
2023-03-04
Earthquakes & Induced quakes- Know the difference
Introduccion
La informacion que usaron para este Kaggle se tomo del USGS(U.S Geological Survey). Este contiene datos e informacion cientifica verificada sobre sucesos geologicos que conciernen a las personas.
El Servicio Geologico Estadounidense es el encargado de monitorear y reportar acerca de terremotos, sus causas, impactos y peligros que estos generan. Tambien envia notificaciones, alertas y demas servicios acerca de terremotos antes que estos ocurran. Se puede encontrar mas informacion en el siguiente link:
Analisis
El autor del Kaggle, cuando extraia la informacion de la pagina, descubrio que no todos los terremotos que se dan se presentan por procesos naturales. Algunos se presentan por actividades humanas como explociones en minas o explociones quimicas. En el periodo de febrero-Marzo de 2019 se dieron mas de 8500 terremotos en el mundo, de los cuales solo el 2% se dieron por actividades en la superrficie.
Al final del trabajo se intento predecir los terremotos y otros tipos de temblores. Asi que el primer paso a seguir es decargar la librerias de informacion.
library(tidyverse) ## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.1.8
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorslibrary(gghighlight)
library(leaflet)
library(IRdisplay)
options(scipen = 999)
options(warn = -1)
library(lubridate)
library(viridis)## Loading required package: viridisLitelibrary(data.table)##
## Attaching package: 'data.table'
##
## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, year
##
## The following objects are masked from 'package:dplyr':
##
## between, first, last
##
## The following object is masked from 'package:purrr':
##
## transposelibrary(caTools)
library(kernlab)##
## Attaching package: 'kernlab'
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## alphalibrary(caret)## Loading required package: lattice
##
## Attaching package: 'caret'
##
## The following object is masked from 'package:purrr':
##
## liftlibrary(ROSE)## Loaded ROSE 0.0-4## Reading in files
list.files("../input")## character(0)earthquake <- read_csv("../all_month.csv")## Rows: 8161 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): magType, net, id, place, type, status, locationSource, magSource
## dbl (12): latitude, longitude, depth, mag, nst, gap, dmin, rms, horizontalE...
## dttm (2): time, updated
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.cols(
.default = col_double(),
time = col_datetime(format = ""),
magType = col_character(),
net = col_character(),
id = col_character(),
updated = col_datetime(format = ""),
place = col_character(),
type = col_character(),
status = col_character(),
locationSource = col_character(),
magSource = col_character()
)## cols(
## .default = col_double(),
## time = col_datetime(format = ""),
## magType = col_character(),
## net = col_character(),
## id = col_character(),
## updated = col_datetime(format = ""),
## place = col_character(),
## type = col_character(),
## status = col_character(),
## locationSource = col_character(),
## magSource = col_character()
## )glimpse(earthquake)## Rows: 8,161
## Columns: 22
## $ time <dttm> 2019-03-11 04:30:18, 2019-03-11 04:24:26, 2019-03-11 …
## $ latitude <dbl> 33.52100, 61.06940, 38.77333, 65.01810, 37.29130, 35.1…
## $ longitude <dbl> -116.7945, -151.8329, -122.7435, -148.7292, -117.5088,…
## $ depth <dbl> 2.82, 105.60, 1.16, 8.60, 6.90, 5.00, 6.20, 48.70, 2.5…
## $ mag <dbl> 0.44, 1.10, 0.73, 1.60, 0.70, 1.80, 1.70, 1.90, 0.56, …
## $ magType <chr> "ml", "ml", "md", "ml", "ml", "mb_lg", "ml", "ml", "md…
## $ nst <dbl> 18, NA, 5, NA, 8, NA, 30, NA, 10, NA, NA, NA, NA, NA, …
## $ gap <dbl> 53.00, NA, 213.00, NA, 176.36, 64.00, 93.14, NA, 77.00…
## $ dmin <dbl> 0.009898, NA, 0.011220, NA, 0.375000, 0.219000, 0.1370…
## $ rms <dbl> 0.1400, 0.2500, 0.0100, 0.5100, 0.2800, 0.5100, 0.3300…
## $ net <chr> "ci", "ak", "nc", "ak", "nn", "us", "nn", "ak", "nc", …
## $ id <chr> "ci38263479", "ak01937u587x", "nc73150441", "ak01937u2…
## $ updated <dttm> 2019-03-11 04:33:55, 2019-03-11 04:27:53, 2019-03-11 …
## $ place <chr> "11km NE of Aguanga, CA", "51km NW of Nikiski, Alaska"…
## $ type <chr> "earthquake", "earthquake", "earthquake", "earthquake"…
## $ horizontalError <dbl> 0.25, NA, 1.55, NA, NA, 1.10, NA, NA, 0.42, 5.80, NA, …
## $ depthError <dbl> 1.02, 1.50, 0.70, 0.20, 24.50, 1.70, 0.70, 1.50, 1.31,…
## $ magError <dbl> 0.129, NA, 0.230, NA, NA, 0.255, NA, NA, NA, 0.080, NA…
## $ magNst <dbl> 11, NA, 2, NA, NA, 4, NA, NA, 1, 15, NA, NA, 11, 60, 1…
## $ status <chr> "automatic", "automatic", "automatic", "automatic", "a…
## $ locationSource <chr> "ci", "ak", "nc", "ak", "nn", "us", "nn", "ak", "nc", …
## $ magSource <chr> "ci", "ak", "nc", "ak", "nn", "us", "nn", "ak", "nc", …head(earthquake)## # A tibble: 6 × 22
## time latitude longi…¹ depth mag magType nst gap dmin
## <dttm> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 2019-03-11 04:30:18 33.5 -117. 2.82 0.44 ml 18 53 0.00990
## 2 2019-03-11 04:24:26 61.1 -152. 106. 1.1 ml NA NA NA
## 3 2019-03-11 04:22:42 38.8 -123. 1.16 0.73 md 5 213 0.0112
## 4 2019-03-11 04:13:42 65.0 -149. 8.6 1.6 ml NA NA NA
## 5 2019-03-11 04:10:13 37.3 -118. 6.9 0.7 ml 8 176. 0.375
## 6 2019-03-11 04:04:12 35.1 -97.6 5 1.8 mb_lg NA 64 0.219
## # … with 13 more variables: rms <dbl>, net <chr>, id <chr>, updated <dttm>,
## # place <chr>, type <chr>, horizontalError <dbl>, depthError <dbl>,
## # magError <dbl>, magNst <dbl>, status <chr>, locationSource <chr>,
## # magSource <chr>, and abbreviated variable name ¹longitudenrow(earthquake)## [1] 8161Resumen sobre las variables y los terminos de los eventos sismicos
Depth- La profundidad en kilometros del evento. Es relativa al nivel del mar (WGS84 geoid) o a la elevacion que presentan las estaciones que toman medidas del terremoto.
depthError- La incertidumbre en la profundidad reportada para el evento.
dmin- La distancia horizontal desde el epicentro hasta la estacion sismologica mas cercana, dada en grados. Un grado es alrededor de 111.2 Km. Entre mas pequeño sea el numero, mas confiable a a ser el calculo de la profundidad.
gap- El espacio mas grande entre estaciones Azimutalmente adyacentes, dado en grados. Entre mas pequeño sea el numero, mas confiable a a ser el calculo de la profundidad. si el espacio sobrepasa los 180 grados, existe mas incertumbre acerca de profundidad y localizacion.
horizontalError- Incertidumbre de la locolizacion reportada, dada en Km
id- Identificador unico del evento.
latitude- Latitud en grados decimales. Valores negativos para latitudes al sur.
locationSource- La fuente que originalmente reporto la localizacion del evento
longitude- Longitud en grados decimales. Valores negativos para longitudes oeste
mag- Magnitud del evento, se mide en base a la fuerza del sismo en su origen
magError- Incertidumbre de la magnitud reportada
magNst- El numero total de estaciones usadas para obtener la magnitud del evento
magSource- La fuente que originalmente reporto la magnitud del evento
magType- El metodo usado para calcular la magnitud del evento
net- La informacion del contribuidor de la informacion del evento
nst- El numero total de estaciones sismicas usadas para determinar la ubicacion del terremoto
place- Descripcion textual de la region geografica mas cercana al sitio donde ocurrio el evento, como una ciudad o un pueblo, en caso que no haya nada cercano a 300 Km, se usa la regionalizacion Flinn-Engdahl (F-E)
rms-El residuo de la media cuadratica del tiempo de viaje, en segundos.
status- Indica si el evento fue revisado por un humano, puede ser automatico o revisado
time- La hora a la que ocurrio el evento, se reporta en milisegundos desde la epoca.
type- Tipo de evento sismico
updated- La hora en que el evento fue actualizado. Se reporta en milisegundos desde la epoca.
Analisis Exploratorio
En kos datos tomados, se presentaro alrededor de 7 sismos con una valoracion mayor a 6 en la escal de Ritcher
options(repr.plot.width=12, repr.plot.height=6)
ggplot(earthquake, aes(time, mag)) +
geom_line(color = 'steelblue')+ theme_bw()+
ggtitle("Earthquake Magnitude")+
xlab('Date')+
ylab('Magnitude')+
theme(plot.title = element_text(hjust = 0.5))
Extraccion de caracteristicas:
Ahora, se extrajeron 2 caracteristicas de los datos:
Localizacion del Sismo (Estado/Pais)
Momento del dia en Horas
De las lineas del siguiente grafico, podemos apreciar el numero maximo de terremotos que ocurrieron a una hora particular.
Se aprecia que medianoche es la hora a la cual el mayor numero de sismos ocurren. Alrededor de las 5 P.M. en cambio, es cuando menos se presentan. Cuando se observa la magnitud maxima por hora, 8-9 A.M. presenta el sismo de mayor impacto, con un 7.
#Location
earthquake$location <- sub('.*,\\s*','', earthquake$place)
#Time of the day in 'Hour'
earthquake$hour <- ymd_hms(earthquake$time)
earthquake$hour <- hour(earthquake$hour)
#Visualizing the number of quakes that have happened at a particular time
earthquake %>%
filter(!is.na(mag))%>%
group_by(hour)%>%
summarise(count = length(id),max_magnitude = max(mag))%>%
ggplot(aes(hour,count, color = max_magnitude))+geom_line()+
scale_color_viridis(direction = -1)+
scale_x_continuous(breaks=seq(0,23,1))+
xlab("Time of the day")+
ylab("Number of earthquakes")+
ggtitle("Number of quakes as per time of the day")+
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
En el mapa que se presenta a continuacion, se aprecian aquellos sitios
donde se presentaron los sismos mas recientes. Ademas, se presentan los
sismos de mayor magnitud mas importantes que se han dado(Entre el 18 y
el 19 de febrero y Marzo)
Place: 115km ESE of Palora, Ecuador Magnitude: 7.5 Time: 2019-02-22 10:17:22
Place: 27km NNE of Azangaro, Peru Magnitude: 7 Time: 2019-03-01 08:50:41
Place: 116km SE of L’Esperance Rock, New Zealand Magnitude: 6.4 Time: 2019-03-06 15:46:14
Place: 49km NW of Namatanai, Papua New Guinea Magnitude: 6.4 Time: 2019-02-17 14:35:55
Place: 28km S of Cliza, Bolivia Magnitude: 6.3 Time: 2019-03-15 05:03:50
Place: 260km SE of Lambasa, Fiji Magnitude: 6.2 Time: 2019-03-10 08:12:25
Place: 140km SSW of Kulumadau, Papua New Guinea Magnitude: 6.1 Time: 2019-03-10 12:48:00
bins=seq(1, 8.0, by=1.0)
palette = colorBin( palette="YlOrBr", domain=earthquake$mag, na.color="transparent", bins=bins)
d=leaflet(earthquake) %>%
addTiles() %>%
setView( lng = 166.45, lat = 21, zoom = 1.25) %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addCircleMarkers(~longitude, ~latitude,
fillColor = ~palette(mag), fillOpacity = 0.7, color="white", radius=3, stroke=FALSE,
popup = paste("Place:", earthquake$place, "<br>",
"Magnitude:", earthquake$mag, "<br>",
"Time:", earthquake$time, "<br>")) %>%
addLegend( pal=palette, values=~mag, opacity=0.9, title = "Magnitude", position = "bottomright" )
htmlwidgets::saveWidget(d, "d.html")
display_html('<iframe src="d.html" width=100% height=450></iframe>')En el Top 30 de las localizaciones donde se han presentado la mayor cantidad se encuentran varios Estados de US asi como paises; aunque alrededor del 35% de los terremotos se presentaron en la costa oeste de los estados unidos y un 27% en Alaska, no superan una magnitud de 1.7, mientras que zonas del hemisferio sur como Chile, presentan menos sismos, pero tienen una magnitud promedio de 4.4
earthquake %>%
group_by(location) %>%
filter(!(is.na(mag)))%>%
summarise(Number_of_quakes = length(location),
Average_Magnitude = mean(mag))%>%
mutate(Percent = round(prop.table(Number_of_quakes)*100,2))%>%
arrange(desc(Number_of_quakes))%>% top_n(25, Number_of_quakes)## # A tibble: 25 × 4
## location Number_of_quakes Average_Magnitude Percent
## <chr> <int> <dbl> <dbl>
## 1 CA 2643 1.02 32.4
## 2 Alaska 2551 1.72 31.3
## 3 Nevada 480 0.776 5.88
## 4 Utah 480 1.39 5.88
## 5 Puerto Rico 269 2.37 3.3
## 6 Hawaii 241 1.83 2.95
## 7 Montana 174 1.27 2.13
## 8 Washington 114 1.06 1.4
## 9 Indonesia 102 4.60 1.25
## 10 California 86 0.980 1.05
## # … with 15 more rowsConocemos que la profundidad es la distancia a la que el terremoto comienza a romper, Asi que a continuacion el autor experimento si a mayor profundidad se podia producir una mayor magnitud. Se creo un grafico de dispersion de profundidad vs magnitud y se encontro que los sismos de menor magnitud suelen ocurrir a mayores profundidades (200 Km o mas), mientras que los de magniudes mayores a 6 suelen estar mas esparcidos.
ggplot(earthquake, aes(mag,depth,color = hour))+
geom_jitter(alpha = 0.5)+
gghighlight(max(depth>200)| max(mag>4))+
scale_x_continuous(breaks=seq(0,9,1))+
scale_color_viridis(direction = -1)+
theme_bw()+
xlab('Magnitude')+
ylab('Depth')+
ggtitle('Magnitude Vs. Depth in Km')+
theme(plot.title = element_text(hjust = 0.5))
Las actividades sismicas tambien pueden presentarse por actividades humanas, como explosiones en canteras, ettc; aun asi, el 98% siguen produciendose por movimientos naturales en la tierra
earthquake %>%
group_by(type)%>%
summarise(count = length(type))%>%
mutate(Percent = prop.table(count*100))%>%
ggplot(aes(type,Percent, fill = type))+
geom_col()+theme_bw()+
xlab('Type of Seismic Activity')+
ylab('Percent')+
ggtitle('Percent of Seismic Activity')+
theme(plot.title = element_text(hjust = 0.5))
Cuando visualizamos las actividades sismicas de acuerdo con su tipo,
podemos ver que gran parte de los movimientos de baja magnitud se
presentan en los Estados Unidos.
pal <- colorFactor(palette = "Accent",domain = earthquake$type)
e=leaflet(earthquake) %>%
addTiles() %>%
setView( lng = -119.417931, lat = 50.778259, zoom = 3) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addCircleMarkers(~longitude, ~latitude,
fillColor = ~pal(type), fillOpacity = 0.5, color="white", radius=3, stroke=FALSE,
popup = paste("Place:", earthquake$place, "<br>",
"Magnitude:", earthquake$mag, "<br>",
"Time:", earthquake$time, "<br>",
"Type:", earthquake$type, "<br>")) %>%
addLegend( pal=pal, values=~type, opacity=0.9, title = "Type", position = "bottomright" )
htmlwidgets::saveWidget(e, "e.html")
display_html('<iframe src="e.html" width=100% height=450></iframe>')Como diferenciar Terremottos de otros movimientos sismicos?
Profundidad
Sabemos que la profundidad es donde el terremoto empieza a moverse. Del grafico de caja que se encuentra abajo, podemos observar que aquellos sismos que no ocurren de forma natural se producen mayormente en la superficie, a diferencia de los naturales, que se producen varios kilometros dentro del planeta.
ggplot(earthquake, aes('',depth, color = type))+geom_boxplot()+
scale_y_continuous(limits = c(0, 100))+theme_bw()+
xlab('Type of Seismic Activity')+
ggtitle('Depth of Seismic Activity')
Magnitud
El siguiente parametro importante es la magnitud. De forma natural, los sismos pueden alcanzar magnitudes entre 0 y 7, mientras los causados por la actividad humana van de 0 a 2.2 en la escala Ritcher.
ggplot(earthquake, aes('',mag, color = type))+geom_boxplot()+
scale_y_continuous(limits = c(0, 8))+theme_bw()+
xlab('Type of Seismic Activity')+
ggtitle('Magnitude of Seismic Activity')
Momento de ocurrencia
Segun se observa, en horas nocturnas como medianoche hasta las 5 A.M. se producen mayormente sismos de actividad humana.
ggplot(earthquake, aes('',hour, color = type))+geom_boxplot()+
scale_y_continuous(limits = c(0, 24))+theme_bw()+
xlab('Type of Seismic Activity')+
ggtitle('Time of Seismic Activity')
Prediccion:
A continuacion, se va a tratar de predecir el tipo de evento sismico por medio de caracteristicas habladas anteriormente.
Primero que nada, se prepara el conjunto de datos para realizar las predicciones. Para los eventos sismicos causados por el hombre, se van a poner todos en una sola categoria ya que tienen un rango de profundidad similar.
#Choosing the relevant columns
quake <- earthquake[,c(2,3,4,5,24,15)]
#Categorising the type
quake$type <- ifelse(quake$type == 'earthquake', 'Earthquake', 'Otherquake')
#quake[,c(1:5)] <- scale(quake[,c(1:5)])
#Converting our target column 'type' into a factor
quake$type <- factor(quake$type)
#Checking and removing any row having NA
colSums(is.na(quake))## latitude longitude depth mag hour type
## 0 0 0 4 0 0quake <- quake[!is.na(quake$mag),]Se usara la maquina de soporte de vectores para hacer predicciones, y se espera un problema de desequilibrio, dado que un 98% de las observaciones son eventos naturales y un 2% artiiciales.
Primero, se entreno al modelo sin tomar en cuenta el desequilibrio, y luego se aplico ROSE library para resolver el problema de equilibrio.
head(quake)## # A tibble: 6 × 6
## latitude longitude depth mag hour type
## <dbl> <dbl> <dbl> <dbl> <int> <fct>
## 1 33.5 -117. 2.82 0.44 4 Earthquake
## 2 61.1 -152. 106. 1.1 4 Earthquake
## 3 38.8 -123. 1.16 0.73 4 Earthquake
## 4 65.0 -149. 8.6 1.6 4 Earthquake
## 5 37.3 -118. 6.9 0.7 4 Earthquake
## 6 35.1 -97.6 5 1.8 4 Earthquaketable(quake$type)##
## Earthquake Otherquake
## 8005 152set.seed(1234)
indices = sample.split(quake$type, SplitRatio = 0.7)
train = quake[indices,]
test = quake[!(indices),]
head(train)## # A tibble: 6 × 6
## latitude longitude depth mag hour type
## <dbl> <dbl> <dbl> <dbl> <int> <fct>
## 1 33.5 -117. 2.82 0.44 4 Earthquake
## 2 61.1 -152. 106. 1.1 4 Earthquake
## 3 38.8 -123. 1.16 0.73 4 Earthquake
## 4 65.0 -149. 8.6 1.6 4 Earthquake
## 5 35.1 -97.6 5 1.8 4 Earthquake
## 6 37.3 -117. 6.2 1.7 3 Earthquaketable(train$type)##
## Earthquake Otherquake
## 5604 106table(test$type)##
## Earthquake Otherquake
## 2401 46set.seed(1234)
#Using RBF Kernel
Model_RBF <- ksvm(type~ ., data = train, scale = FALSE, kernel = "rbfdot")
Eval_RBF <- predict(Model_RBF, test[,-6])
#confusion matrix - RBF Kernel
confusionMatrix(Eval_RBF,test$type)## Confusion Matrix and Statistics
##
## Reference
## Prediction Earthquake Otherquake
## Earthquake 2401 46
## Otherquake 0 0
##
## Accuracy : 0.9812
## 95% CI : (0.975, 0.9862)
## No Information Rate : 0.9812
## P-Value [Acc > NIR] : 0.5391
##
## Kappa : 0
##
## Mcnemar's Test P-Value : 0.00000000003247
##
## Sensitivity : 1.0000
## Specificity : 0.0000
## Pos Pred Value : 0.9812
## Neg Pred Value : NaN
## Prevalence : 0.9812
## Detection Rate : 0.9812
## Detection Prevalence : 1.0000
## Balanced Accuracy : 0.5000
##
## 'Positive' Class : Earthquake
## En la prediccion de arriba, se observa con claridad el problema. Todos los terremotos naturales fueron predecidos adecuadamente, mientras que los artificiales no pudieron ser predecidos en lo mas minimo.
Solucionando el desbalanceo
Se va a aplicar el paquete de ROSE (Randomly Over Sampling Examples) para solucionar el problema con el balanceo mediate sobremuestreo y submuestreo.
set.seed(1234)
over <- ovun.sample(type~., data = train, method="over", N= 11640)$dataset.seed(1234)
#Using RBF Kernel
Model_RBF <- ksvm(type~ ., data = over, scale = FALSE, kernel = "rbfdot")
Eval_RBF <- predict(Model_RBF, test[,-6])
#confusion matrix - RBF Kernel
confusionMatrix(Eval_RBF,test$type)## Confusion Matrix and Statistics
##
## Reference
## Prediction Earthquake Otherquake
## Earthquake 2285 4
## Otherquake 116 42
##
## Accuracy : 0.951
## 95% CI : (0.9416, 0.9592)
## No Information Rate : 0.9812
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.3941
##
## Mcnemar's Test P-Value : <0.0000000000000002
##
## Sensitivity : 0.9517
## Specificity : 0.9130
## Pos Pred Value : 0.9983
## Neg Pred Value : 0.2658
## Prevalence : 0.9812
## Detection Rate : 0.9338
## Detection Prevalence : 0.9354
## Balanced Accuracy : 0.9324
##
## 'Positive' Class : Earthquake
## Ahora se puede obsevar un aumento en la especifidad, de alrededor del 85%.