R Visualization Project

Exploration of Earthquakes Recorded in 2023

In 2023, a large range of Earthquakes were recorded. In this dataset, the location (long/lat and descriptive), time, magnitude, and more was entered. Through the next visualizations, I will explore the quake types, the locations, the magnitudes, and the season in which the quake occurred.

setwd("C://Users//ashkl//OneDrive//Documents//Loyola//data visualization")
df<- read.csv("earthquakes_2023_global.csv" , na.strings=c(NA, ""))
library(dplyr)
library(ggplot2)
library(leaflet)
library(maps)
library(lubridate)
library(scales)
library(plotly)

Basic Statistics

Navigate through the below tabs to view some basic statistics of the data.

Summary

summary(df)

##      time              latitude         longitude           depth       
##  Length:26642       Min.   :-65.850   Min.   :-180.00   Min.   : -3.37  
##  Class :character   1st Qu.: -6.415   1st Qu.:-149.61   1st Qu.: 10.00  
##  Mode  :character   Median : 18.884   Median : -64.81   Median : 22.00  
##                     Mean   : 16.853   Mean   : -11.49   Mean   : 67.49  
##                     3rd Qu.: 41.828   3rd Qu.: 126.97   3rd Qu.: 66.83  
##                     Max.   : 86.594   Max.   : 180.00   Max.   :681.24  
##                                                                         
##       mag          magType               nst              gap       
##  Min.   :2.600   Length:26642       Min.   :  0.00   Min.   :  8.0  
##  1st Qu.:3.220   Class :character   1st Qu.: 19.00   1st Qu.: 73.0  
##  Median :4.300   Mode  :character   Median : 30.00   Median :111.0  
##  Mean   :4.007                      Mean   : 42.57   Mean   :124.9  
##  3rd Qu.:4.500                      3rd Qu.: 52.00   3rd Qu.:165.0  
##  Max.   :7.800                      Max.   :423.00   Max.   :350.0  
##                                     NA's   :1415     NA's   :1417   
##       dmin             rms             net                 id           
##  Min.   : 0.000   Min.   :0.0100   Length:26642       Length:26642      
##  1st Qu.: 0.612   1st Qu.:0.4100   Class :character   Class :character  
##  Median : 1.579   Median :0.5900   Mode  :character   Mode  :character  
##  Mean   : 2.693   Mean   :0.5816                                        
##  3rd Qu.: 3.172   3rd Qu.:0.7500                                        
##  Max.   :50.820   Max.   :1.8800                                        
##  NA's   :1866                                                           
##    updated             place               type           horizontalError 
##  Length:26642       Length:26642       Length:26642       Min.   : 0.000  
##  Class :character   Class :character   Class :character   1st Qu.: 4.140  
##  Mode  :character   Mode  :character   Mode  :character   Median : 7.060  
##                                                           Mean   : 7.017  
##                                                           3rd Qu.: 9.730  
##                                                           Max.   :99.000  
##                                                           NA's   :1549    
##    depthError        magError          magNst          status         
##  Min.   : 0.000   Min.   :0.0000   Min.   :  0.00   Length:26642      
##  1st Qu.: 1.848   1st Qu.:0.0800   1st Qu.: 10.00   Class :character  
##  Median : 2.019   Median :0.1110   Median : 18.00   Mode  :character  
##  Mean   : 4.475   Mean   :0.1227   Mean   : 33.32                     
##  3rd Qu.: 6.669   3rd Qu.:0.1500   3rd Qu.: 36.00                     
##  Max.   :60.670   Max.   :4.4900   Max.   :884.00                     
##                   NA's   :1672     NA's   :1577                       
##  locationSource      magSource        
##  Length:26642       Length:26642      
##  Class :character   Class :character  
##  Mode  :character   Mode  :character  
##                                       
##                                       
##                                       
## 

Structure

str(df)

## 'data.frame':    26642 obs. of  22 variables:
##  $ time           : chr  "2023-01-01T00:49:25.294Z" "2023-01-01T01:41:43.755Z" "2023-01-01T03:29:31.070Z" "2023-01-01T04:09:32.814Z" ...
##  $ latitude       : num  52.1 7.14 19.16 -4.78 53.4 ...
##  $ longitude      : num  178.5 126.7 -66.5 102.8 -166.9 ...
##  $ depth          : num  82.8 79.2 24 63.8 10 ...
##  $ mag            : num  3.1 4.5 3.93 4.3 3 2.8 4.1 4.1 2.7 4.1 ...
##  $ magType        : chr  "ml" "mb" "md" "mb" ...
##  $ nst            : num  14 32 23 17 19 19 15 40 30 27 ...
##  $ gap            : num  139 104 246 187 190 127 87 81 154 119 ...
##  $ dmin           : num  0.87 1.152 0.848 0.457 0.4 ...
##  $ rms            : num  0.18 0.47 0.22 0.51 0.31 0.18 0.15 0.32 0.14 0.45 ...
##  $ net            : chr  "us" "us" "pr" "us" ...
##  $ id             : chr  "us7000j5a1" "us7000j3xk" "pr2023001000" "us7000j3xm" ...
##  $ updated        : chr  "2023-03-11T22:51:52.040Z" "2023-03-11T22:51:45.040Z" "2023-03-11T22:51:29.040Z" "2023-03-11T22:51:45.040Z" ...
##  $ place          : chr  "Rat Islands, Aleutian Islands, Alaska" "23 km ESE of Manay, Philippines" "Puerto Rico region" "99 km SSW of Pagar Alam, Indonesia" ...
##  $ type           : chr  "earthquake" "earthquake" "earthquake" "earthquake" ...
##  $ horizontalError: num  8.46 5.51 0.91 10.25 1.41 ...
##  $ depthError     : num  21.21 7.45 15.95 6.58 2 ...
##  $ magError       : num  0.097 0.083 0.09 0.238 0.085 0.06 0.213 0.095 0.112 0.123 ...
##  $ magNst         : num  14 43 16 5 18 36 6 34 35 18 ...
##  $ status         : chr  "reviewed" "reviewed" "reviewed" "reviewed" ...
##  $ locationSource : chr  "us" "us" "pr" "us" ...
##  $ magSource      : chr  "us" "us" "pr" "us" ...

Columns With NA

colSums(is.na(df))

##            time        latitude       longitude           depth             mag 
##               0               0               0               0               0 
##         magType             nst             gap            dmin             rms 
##               0            1415            1417            1866               0 
##             net              id         updated           place            type 
##               0               0               0            1608               0 
## horizontalError      depthError        magError          magNst          status 
##            1549               0            1672            1577               0 
##  locationSource       magSource 
##               0               0

Count of Quake Causes

eq_cause_count <- df %>%
  select(type) %>%
  group_by(type) %>%
  tally() %>%
  data.frame()
eq_cause_count

##                type     n
## 1         Landslide     2
## 2        earthquake 26428
## 3         explosion     3
## 4         ice quake    26
## 5         landslide     1
## 6  mining explosion   167
## 7      quarry blast     2
## 8 volcanic eruption    13

#initial setup of df
df<- na.omit(df)
df <- subset(df, select = -c(status, updated))
df$country <- gsub(".*, ", "", df$place)

#create df for count of quake types to be used in bar chart
eq_cause_count <- df %>%
  select(type) %>%
  group_by(type) %>%
  tally() %>%
  data.frame()

eq_cause_other_pct<- eq_cause_count[-c(1), ]
eq_cause_other_pct$pct <- paste0(round(eq_cause_other_pct$n/sum(eq_cause_other_pct$n) *100, digits =1),"%")

#create dfs for only mining and quarry blast quake types and for volcano types to be used for map
eq_mining <- df[df$type %in% c('mining explosion', 'quarry blast'),]
eq_volcano<- df[df$type %in% c('volcanic eruption'),]

#create df for only earthquakes and one for only non-earth type quakes to be used in dot plot map of all quakes
df$type <- as.factor(df$type)
eq_other<- df[!df$type %in% c('earthquake'),]
earthquake_only <-  df[df$type %in% c('earthquake'),]
world_map <- map_data("world")

#create dfs to categorize quake month into seasons with respect to longitude because seasons differ by north/south hemispheres and by equator location then bind them back together
df$time <- as.Date(df$time)
df$season <- month(df$time)

equator_quakes <- df[df$longitude <= 25 & df$longitude >= -25,]
equator_quakes$season <- "Summer"

NHem_quakes <- df[df$longitude > 25,]
NHem_quakes["season"][NHem_quakes["season"] == 1 | NHem_quakes["season"] == 2 | NHem_quakes["season"] == 12] <- "Winter"
NHem_quakes["season"][NHem_quakes["season"] == 3 | NHem_quakes["season"] == 4 | NHem_quakes["season"] == 5] <- "Spring"
NHem_quakes["season"][NHem_quakes["season"] == 6 | NHem_quakes["season"] == 7 | NHem_quakes["season"] == 8] <- "Summer"
NHem_quakes["season"][NHem_quakes["season"] == 9 | NHem_quakes["season"] == 10 | NHem_quakes["season"] == 11] <- "Fall"
SHem_quakes <- df[df$longitude < -25,]
SHem_quakes["season"][SHem_quakes["season"] == 1 | SHem_quakes["season"] == 2 | SHem_quakes["season"] == 12] <- "Summer"
SHem_quakes["season"][SHem_quakes["season"] == 3 | SHem_quakes["season"] == 4 | SHem_quakes["season"] == 5] <- "Fall"
SHem_quakes["season"][SHem_quakes["season"] == 6 | SHem_quakes["season"] == 7 | SHem_quakes["season"] == 8] <- "Winter"
SHem_quakes["season"][SHem_quakes["season"] == 9 | SHem_quakes["season"] == 10 | SHem_quakes["season"] == 11] <- "Spring"

quakes_w_season <- do.call("rbind", list(equator_quakes, NHem_quakes, SHem_quakes)) 

#get the top 10 countries/states with quakes
place_count <- quakes_w_season %>%
  select(season, country) %>%
  group_by(country) %>%
  tally() %>%
  data.frame()

place_count <- place_count[order(place_count$n, decreasing = TRUE), ]
place_count <- place_count[1:10, ]

#get seasons counts of only the top 10
keep<- place_count$country
place_count_season <- subset(quakes_w_season, country %in% keep)
place_count_season <- place_count_season %>%
  select(season, country) %>%
  group_by(country, season) %>%
  tally() %>%
  data.frame()

#create season levels for heatmap
mylevels<- c("Spring", "Summer", "Fall", "Winter")
place_count_season$season <- factor(place_count_season$season, levels = mylevels)

#create df for the magnitudes of places in top 10, then factor them
df_magnitudes <- subset(df, country %in% keep)

df_magnitudes <- df_magnitudes %>%
  select(country, mag) %>%
  mutate(mag_levels = ifelse(mag <= 6, "Minor", ifelse(mag > 6 & mag < 7, "Moderate", ifelse(mag >= 7 & mag < 8, "Major", "Massive")))) %>%
  group_by(country, mag_levels) %>%
  tally() %>%
  group_by(country) %>%
  mutate (pct_of_total = round(100*n/sum(n), 1)) %>%
  ungroup() %>%
  data.frame()

df_magnitudes$mag_levels = factor(df_magnitudes$mag_levels, levels = c("Minor", "Moderate", "Major", "Massive"))

Bar Chart

Due to the vast majority of quakes being naturally-occurring earthquakes, included in the below are only the non “Eath” type quakes. As you can see, the majority of these quakes were caused by mining explosions and volcanoes. Were these located in one area?

ggplot(eq_cause_other_pct, aes(x=reorder(type,-n), y=n)) + 
  geom_histogram(colour = "black", fill="brown", stat="identity") +
  labs(title = "Non-normal Earthquakes in 2023", x="Quake Cause", 
       y="Quake Count") +   
  geom_text(aes(label=eq_cause_other_pct$pct), vjust=-0.3, size=3.5) +
  theme_light() +
  theme(plot.title = element_text(hjust=.5, size = 16))

Interactive Maps

To answer the question of where the volcano and mining quakes happened, navigate through the below tabs to view interactive maps of where each type occurred.

Map of Volcanoes

As we can see if we zoom out just a little, all of the volcano quakes were recorded in ocean just North of Northern Mariana Islands.

#Volcanoes
icon.glyphicon <- makeAwesomeIcon(icon = 'fire', markerColor = 'red', iconColor = 'white')

map_volcanoes <- leaflet() %>%
  addTiles() %>%
  addAwesomeMarkers(lng = eq_volcano$longitude, lat = eq_volcano$latitude,
                    icon = icon.glyphicon,
                    popup = paste("Magnitude:", eq_volcano$mag),
                    label = eq_volcano$place)
map_volcanoes

Map of Mining/Quarry Blasts

Surprising, but perhaps not so surprising that all (minus one exception) quakes caused by mining and quarry blasts recorded were located in the United States. This could be due to other countries not classifying the cause, but it could also be that the explosions in other countries were not strong enough to be picked up by a local seismograph.

#quarry/mining blasts
icon.glyphicon <- makeAwesomeIcon(icon = 'star', markerColor = 'red', iconColor = 'white')

map_mining <- leaflet() %>%
  addTiles() %>%
  addAwesomeMarkers(lng = eq_mining$longitude, lat = eq_mining$latitude,
                    icon = icon.glyphicon,
                    popup = paste(eq_mining$place, "| Magnitude:", eq_mining$mag),
                    label = eq_mining$type
                    )
map_mining

Dot Plot

A dot plot was used to show the location of every quake recorded in 2023. Why a dot plot? Because my poor computer could not load an interactive leaftlet plot and because with a dot plot, I can show where each type of quake occurred.

map_dot<- ggplot() +
  coord_fixed() +
  xlab("") +
  ylab("")+
  geom_polygon(data=world_map, aes(x = long, y = lat, group = group),
                 colour = 'light grey', fill = 'light grey') +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white", colour = "white"),
        axis.line = element_line(colour = "white"), 
        axis.ticks = element_blank(), axis.text.x = element_blank(),
        axis.text.y = element_blank()) +  
  geom_point(data = earthquake_only, aes(x = longitude, y = latitude), shape = 21, size = 2, fill = 'white', color = 'black') +
  geom_point(data = eq_other, aes(x = longitude, y = latitude, color = type, shape = type, fill = type), size = 4) +
  scale_shape_manual(values=c(8, 18, 9, 15, 2)) +
  labs(title = "Dot Plot with Map of Quakes and Types 2023", x="Latitude", y = "Longitude") +
  theme(plot.title = element_text(hjust = .5))

map_dot

Heatmap

This was an interesting find. I created a new dataframe and column for each season (respective of the logitude). The Northern Hemisphere seasons were applied to longitude > 25. Equater season (summer) was applied to longitude <= 25 through longitude >= -25. The Southern Hemisphere seasons were applide to longitude < -25. I found that for most coutries/states, quakes were more common in Winter. Now, in the below heatmap of the top 10 places where quakes occurred, this was not always true, but this was the general trend.

map_heat<- ggplot(place_count_season, aes(x = season, y = country, fill = n))+
  geom_tile(color = "black") + 
  geom_text(aes(label = n)) + 
  coord_equal(ratio = 1) +
  labs(title = "Heatmap: Count of Quakes by Season by Country/State", x = "Season Quake Occurred", y = "State or Country of Quake", fill = "Quake Count") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_x_discrete(limits = rev(levels(place_count_season$season))) +
  scale_fill_continuous(low = "white", high = "blue") +
  guides(fill = guide_legend(reverse = TRUE, override.aes = list(colour = "black"))) 


ggplotly(map_heat, tooltip = c("n", "country", "season")) %>%
  style(hoverlabel = list(bgcolor = "white"))

Donut Charts for Top 10 Quakes

For the top 10 countries where quakes occurred, donut charts were created for the magnitude strength. I have separated the top 10 into separate charts by region. Navigate through the below tabs to view the magnitudes recorded in the United States, Oceania, Asia, and Europe.Magnitude strength is listed as “Minor,” “Moderate,” and “Major.”

United States Magnitudes

US_States <- plot_ly(hole = .7) %>%
  layout(title = "Magnitude of Quakes in the United States") %>%
  add_trace(data = df_magnitudes[df_magnitudes$country == "Alaska",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Alaska", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Alaska<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>") %>%
add_trace(data = df_magnitudes[df_magnitudes$country == "CA",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "CA", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: California<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>",
domain = list(x = c(0.16, 0.84), y = c(0.16, 0.84))) %>%
add_trace(data = df_magnitudes[df_magnitudes$country == "Puerto Rico",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Puerto Rico", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Puerto Rico<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>",
          domain = list(x = c(0.27, 0.73), y = c(0.27, 0.73)))
US_States

Oceania Magnitudes

Oceania <- plot_ly(hole = .7) %>%
  layout(title = "Magnitude of Quakes in the Oceania") %>%
  add_trace(data = df_magnitudes[df_magnitudes$country == "Tonga",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Tonga", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Tonga<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>") %>%
add_trace(data = df_magnitudes[df_magnitudes$country == "Papua New Guinea",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Papua New Guinea", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Papua New Guinea<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>",
domain = list(x = c(0.16, 0.84), y = c(0.16, 0.84))) 
Oceania

Asia Magnitudes

Asia <- plot_ly(hole = .7) %>%
  layout(title = "Magnitude of Quakes in Asia") %>%
  add_trace(data = df_magnitudes[df_magnitudes$country == "Indonesia",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Indonesia", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Indonesia<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>") %>%
add_trace(data = df_magnitudes[df_magnitudes$country == "Japan",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Japan", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Japan<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>",
domain = list(x = c(0.16, 0.84), y = c(0.16, 0.84))) %>%
add_trace(data = df_magnitudes[df_magnitudes$country == "Japan region",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Japan region", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Japan Region<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>",
          domain = list(x = c(0.27, 0.73), y = c(0.27, 0.73))) %>%
add_trace(data = df_magnitudes[df_magnitudes$country == "Philippines",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Philippines", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Philippines<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>",
          domain = list(x = c(0.345, 0.655), y = c(0.345, 0.655)))
Asia

Europe Magnitudes

Europe <- plot_ly(hole = .7) %>%
  layout(title = "Magnitude of Quakes in Europe") %>%
  add_trace(data = df_magnitudes[df_magnitudes$country == "Turkey",], 
            labels = ~mag_levels, 
            values = ~df_magnitudes[df_magnitudes$country == "Turkey", "n"], 
            type = "pie", 
            textposition = "inside", 
            hovertemplate = "Location: Turkey<br>Magnitude: %{label}<br>Percent: %{percent}<br>Quake Count: %{value}<extra></extra>")
Europe

Conclusions

Earthquakes happen all over the world every day. While more uncommon, other quake types such as mining explosions, volcanoes, and ice quakes occur. Most, if not all, natural earthquakes happen around the earth’s technonic plates. As a general rule, more earthquakes occur in the cold winter months; this could be due the natural expansion and contraction of heating and cooling, but more research should be done to confirm this. The vast majority of earthquakes that occurred in 2023 were minor quakes with a magnitude <= 6 on the Richter scale.