INTRODUCCION

Los datos son relevantes para evaluar el impacto y riesgos sísmicos, llevan a cabo investigaciones sobre las causas y efectos de un terremoto. Se usaron los datos de (USGS) U.S Geology survey.

Este estudio surge por el hecho de que todos los terremotos no son causados por la naturaleza, ya que se logró evidenciar que alrededor del 2% de los terremotos mensuales en el 2019 fueron causador por (explosiones de canteras y actividades relacionadas con explosiones)

IMPORTAR LIBRERIAS

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 errors
library(gghighlight)
library(leaflet)
library(IRdisplay)
options(scipen = 999)
options(warn = -1)
library(lubridate)
library(viridis)
## Loading required package: viridisLite
library(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':
## 
##     transpose
library(caTools)
library(kernlab)
## 
## Attaching package: 'kernlab'
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     alpha
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(ROSE)
## Loaded ROSE 0.0-4
list.files(path = "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.
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", …

GRAFICO 1

DATOS

En la primera tabla se aprecia la organización de la recolección de los datos en la que tendremos datos como: la profundidad del evento en kilómetros,incertidumbre en la profundidad reportada ,distancia horizontal desde el epicentro ,incertidumbre de la localización informada ,identificador único del suceso, grados decimales de latitud ,magnitud del sismo y su incertidumbre ,número de estaciones sísmicas usadas para calcular la magnitud etc…

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 ¹​longitude

GRAFICO 2

EARTHQUAKE MAGNITUDE

en esta grafica se puede observar la relación (día- magnitud) es decir la magnitud de los sismos entre febrero y marzo del 2019, se logra evidenciar que alrededor del 18 y 25 de mayo gubo un sismo de magnitud mayor a 6 .

options(repr.plot.width=12, repr.plot.height=6)

ggplot(earthquake, aes(time, mag)) + 
geom_line(color = 'green')+ theme_bw()+ 
ggtitle("Earthquake Magnitude")+
xlab('Date')+
ylab('Magnitude')+
theme(plot.title = element_text(hjust = 0.5))

GRAFICO 3

NÚMERO DE TERREMOTOS SEGÚN LA HORA DEL DÍA

en este grafico de líneas se observan el número de terremotos que se presentan en ciertas horas del día. De lo que se puede deducir que a la media noche se producen la mayor cantidad de terremotos y las horas en las que menos hay terremotos son entre las 5 y las 7 de la noche. Y alrededor de las 8 y las 9 se evidencia el terremoto de mayor magnitud

#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))

GRAFICO 4

En este grafico se puede ver el numero de terremotos en cada pais, ademas de una media de la magnitud de los terremotos en cada pais, ademas de los porcentajes a nivel mundial…

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 rows

GRAFICO 5

MAGNITUD Y PROFUNDIDAD EN KM

Grafico de dispersión magnitud vs profundidad, para saber si la magnitud es directamente proporcional con la profundidad. Y se logró evidenciar que los terremotos de menor magnitud suelen darse 200km de profundidad. Y un numero contado de sismos de menor magnitud tienen epicentro en
altas profundidades.

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))

GRAFICO 6

PORCENTAJES DE ACTIVIDAD SISMICA

En el diagrama de caja vemos que la mayoría de las actividades sísmicas no relacionadas con la naturaleza se producen sobre todo en la superficie de la Tierra. Por otro lado, los terremotos se producen sobre todo entre 10 y 20 km por debajo de la superficie y pueden llegar a producirse hasta varios kilómetros por debajo de la superficie de la Tierra. Pero esto no suele ocurrir con otro tipo de actividades sísmicas

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))

GRAFICO 7

PROFUNDIDAD DE ACTIVIDAD SISMICA

en este diagrama de caja se logra ver que los sismos natural naturales puede tener altas profundidades y hay gran cantidad de datos atipicos, a comparacion de otros eventos que solo ocurren a profundidades minimas.

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')

## GRAFICO 7

MAGNITUD DE ACTIVIDAD SISMICA

en este diagrama de caja vemos la relevancia de la magnitud de un terremoto natural puede tener una magnitud de entre 0 y 7 grados. La magnitud de los terremotos provocados por el hombre es mucho menor y suele oscilar entre 0 y 2,2 en la escala de Richter.

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')

GRAFICO 8

TIEMPO DE ACTIVIDAD SISMICA

En etse grafico de caja se observa que los terremotos en determinadas horas pueden deberse a una explosion quimica de la zona ya que segun las estadisticas de media noche a 5 AM son pocas los registros de sismos naturales.

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')

GRAFICO 9

PREDICCION

En lo siguiente se intento predecir que habia causado el terremoto teniendo en cuenta factores como la profundidad, maginitud y ubicacion ya que la actividad sismica causada naturalmente y provocada, tiene rangos estimados para los factores mencionados

#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         0
quake <- quake[!is.na(quake$mag),]
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 Earthquake
table(quake$type)
## 
## Earthquake Otherquake 
##       8005        152
# splitting the data between train and test
set.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 Earthquake
table(train$type)
## 
## Earthquake Otherquake 
##       5604        106
table(test$type)
## 
## Earthquake Otherquake 
##       2401         46

SIN TRATAR EL DESEQUILIBRIO DE CLASES

En esta prediccion podemos ver claramente el fallo. Todos los terremotos se predicen correctamente, lo que hace que la sensibilidad (verdaderos positivos) sea muy alta. Sin embargo, ninguno de los “otros terremotos” se predice correctamente, lo que lleva a una Especificidad (Verdaderos Negativos) de 0.

set.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      
##

TRATANDO EL DESEQUILIBRIO DE CLASES

Se puede ver claramente un aumento significativo de la Especificidad. 47 de 55 Otros terremotos en nuestros datos de prueba se predicen correctamente aunque 150 terremotos se predicen incorrectamente. La introducción del sobremuestreo ha ayudado mucho a tratar el desequilibrio de clases y ha aumentado significativamente nuestra Especificidad hasta el 85%.

set.seed(1234)
over <- ovun.sample(type~., data = train, method="over", N= 11640)$data
set.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         
##

MAPA DE LOS SISMOS CON SUS MAGNITUDES EN TODO EL MUNDO

En este mapa se logra observar que muchos de los terremtos que suceden en america ocurren en el cinturon de fuego del pacifico y a juzgar por la escala de colores que indican la magnitud del sismo tiene un rango de magnitud de ente (4-5)

#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>')

Mapa 1

MAPA DE UBICACION DE SISMOS TENIENDO EN CUENTA QUE LOS PRODUJO

Teniendo en cuenta que los sismos pueden ser inducidos o pueden ser generado por tectonica de placas, en este mapa se logro evidenciar por colores algunos de los sismos que se han dado por explosiones de canteras y que los sismos que han sido ocacionados por tectonica son mayoria.

#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>')

Mapa 2