Introduction

Earthquakes around the world pose an increasing risk hazard given the growing population in high risk areas. Terrains with beautiful vistas in first world countries may be at risk of mud slides during earthquake shaking events. Populations that live along the sea are at risk due to the aftermath of tsunami’s after a major earthquake event. Farming communities that live in water rich valley floors, frequently along fault zones may be at risk. Crowded cities in overpopulated third world countries are most at risk due to outdated or worse no building codes in overcrowded living conditions with poor to no access to health care coverage. There were 6 quakes greater than magnitude 5 in the past 30 days along the west coast of South America Recent quakes in Chile.

The development of this widget is presented in the sequential steps taken to select the correct base map, to select the correct map view, to select the correct colors to convey the story effectively. The analysis takes the reader through the tough process of why certain steps were taken to build an interactive world earthquakes map in leaflet.

Image Credit: istock.com/Petrovich.

Data Preparation

Base Map Preparation for the interactive earthquake app. Selecting the most appropriate base map to match the purpose of the study is an important part of creating a leaflet interactive map. I decided on Esri’s World Image as a base as it has real-world colors, has clearly defined fault lines in the ocean that continue onto land, where some of the weaker earthquakes were located and it references the reader to actual continents which makes seeing the earthquakes nearby more of a concern.

#* Load the required packages*

library(leaflet)
library(here)
library(readr)
library(tidyverse)
library(htmlwidgets)

set.seed(122) #* to ensure you get the same result when you run the code in R*

To access the data for the interactive map creation you must read in the earthquake data and create a new map object with magnitude values of interest, in this case magnitude greater than 4.5 magnitude. I decided to focus on 50 earthquakes greater than 4.5 magnitude as high as 6.0 magnitude.

#### *Read in the earthquakes greater than 4.5 magnitude data and view the head*

quakes_greater4.5 <- read_csv(here("data","4.5_month.csv"))
head(quakes_greater4.5)
## # A tibble: 6 × 22
##   time                latitude longitude depth   mag magType   nst   gap  dmin
##   <dttm>                 <dbl>     <dbl> <dbl> <dbl> <chr>   <dbl> <dbl> <dbl>
## 1 2023-06-09 22:46:20    16.5      146.  145.    4.5 mb         51    80  1.21
## 2 2023-06-09 22:25:00    -2.81     129.   10     4.7 mb         43    43  1.83
## 3 2023-06-09 21:33:05   -16.7     -174.   92.8   4.8 mb         49    49  4.08
## 4 2023-06-09 21:21:42   -60.5      160.   10     5.9 mww        63    56  6.06
## 5 2023-06-09 21:04:10   -23.3      171.   10     4.8 mb         28    52  3.02
## 6 2023-06-09 20:43:28   -47.3      -12.7  10     5.5 mb        111    60 16.5 
## # ℹ 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>
#* Sample the data set for 50 observations*
quakes1=quakes_greater4.5[sample(nrow(quakes_greater4.5),50),]
str(quakes1) #* Note there are now only 50 observations and 22 variables*
## tibble [50 × 22] (S3: tbl_df/tbl/data.frame)
##  $ time           : POSIXct[1:50], format: "2023-06-07 15:51:24" "2023-05-27 16:49:22" ...
##  $ latitude       : num [1:50] -56.092 60.277 0.966 -23.048 -6.596 ...
##  $ longitude      : num [1:50] -27.3 -153 126.6 170.2 147.1 ...
##  $ depth          : num [1:50] 35 138.7 59.5 10 52.1 ...
##  $ mag            : num [1:50] 4.9 4.6 4.5 5.2 4.9 4.5 5 5.1 4.8 4.6 ...
##  $ magType        : chr [1:50] "mb" "mb" "mb" "mb" ...
##  $ nst            : num [1:50] 30 91 25 52 83 24 82 41 85 27 ...
##  $ gap            : num [1:50] 124 38 77 76 73 102 66 48 37 122 ...
##  $ dmin           : num [1:50] 5.565 0.18 0.781 2.577 2.791 ...
##  $ rms            : num [1:50] 0.49 0.52 0.58 0.96 0.68 0.7 0.76 0.74 0.7 0.78 ...
##  $ net            : chr [1:50] "us" "us" "us" "us" ...
##  $ id             : chr [1:50] "us7000k6xd" "us7000k4hp" "us7000k6dz" "us7000k6qq" ...
##  $ updated        : POSIXct[1:50], format: "2023-06-07 16:25:56" "2023-06-05 17:54:18" ...
##  $ place          : chr [1:50] "South Sandwich Islands region" "73 km E of Port Alsworth, Alaska" "Molucca Sea" "280 km E of Vao, New Caledonia" ...
##  $ type           : chr [1:50] "earthquake" "earthquake" "earthquake" "earthquake" ...
##  $ horizontalError: num [1:50] 12.53 4.23 8.32 6.64 10 ...
##  $ depthError     : num [1:50] 1.92 4.79 8.82 1.86 7.16 ...
##  $ magError       : num [1:50] 0.125 0.137 0.15 0.092 0.051 0.122 0.053 0.114 0.022 0.136 ...
##  $ magNst         : num [1:50] 21 16 13 39 123 20 113 25 642 17 ...
##  $ status         : chr [1:50] "reviewed" "reviewed" "reviewed" "reviewed" ...
##  $ locationSource : chr [1:50] "us" "us" "us" "us" ...
##  $ magSource      : chr [1:50] "us" "us" "us" "us" ...

Analysis

Next I created three earthquake magnitude types by ranges including weak (4-5 magnitude), moderate(5-6 magnitude), strong(6-7 magnitude). I created an intuitive name for quakes by map range then took a look using the str() function to ensure the code worked correctly and a new column called maprange was successfully created.

#* Create an earthquake magnitude range to define the quake strength*
#* Create an object mag range; add variables quake light, quake moderate or quake strong
library(readr)

quakes1$magrange = cut(quakes1$mag, 
                                 breaks = c(4, 5, 6, 7), right = FALSE, #* to exclude the right hand mag ranges*
                                 labels = c("Weak[4-5)","Moderate[5-6)", "Strong[6-7)"))
str(quakes1) #* To view new column magrange with values corresponding to mag ranges defined*
## tibble [50 × 23] (S3: tbl_df/tbl/data.frame)
##  $ time           : POSIXct[1:50], format: "2023-06-07 15:51:24" "2023-05-27 16:49:22" ...
##  $ latitude       : num [1:50] -56.092 60.277 0.966 -23.048 -6.596 ...
##  $ longitude      : num [1:50] -27.3 -153 126.6 170.2 147.1 ...
##  $ depth          : num [1:50] 35 138.7 59.5 10 52.1 ...
##  $ mag            : num [1:50] 4.9 4.6 4.5 5.2 4.9 4.5 5 5.1 4.8 4.6 ...
##  $ magType        : chr [1:50] "mb" "mb" "mb" "mb" ...
##  $ nst            : num [1:50] 30 91 25 52 83 24 82 41 85 27 ...
##  $ gap            : num [1:50] 124 38 77 76 73 102 66 48 37 122 ...
##  $ dmin           : num [1:50] 5.565 0.18 0.781 2.577 2.791 ...
##  $ rms            : num [1:50] 0.49 0.52 0.58 0.96 0.68 0.7 0.76 0.74 0.7 0.78 ...
##  $ net            : chr [1:50] "us" "us" "us" "us" ...
##  $ id             : chr [1:50] "us7000k6xd" "us7000k4hp" "us7000k6dz" "us7000k6qq" ...
##  $ updated        : POSIXct[1:50], format: "2023-06-07 16:25:56" "2023-06-05 17:54:18" ...
##  $ place          : chr [1:50] "South Sandwich Islands region" "73 km E of Port Alsworth, Alaska" "Molucca Sea" "280 km E of Vao, New Caledonia" ...
##  $ type           : chr [1:50] "earthquake" "earthquake" "earthquake" "earthquake" ...
##  $ horizontalError: num [1:50] 12.53 4.23 8.32 6.64 10 ...
##  $ depthError     : num [1:50] 1.92 4.79 8.82 1.86 7.16 ...
##  $ magError       : num [1:50] 0.125 0.137 0.15 0.092 0.051 0.122 0.053 0.114 0.022 0.136 ...
##  $ magNst         : num [1:50] 21 16 13 39 123 20 113 25 642 17 ...
##  $ status         : chr [1:50] "reviewed" "reviewed" "reviewed" "reviewed" ...
##  $ locationSource : chr [1:50] "us" "us" "us" "us" ...
##  $ magSource      : chr [1:50] "us" "us" "us" "us" ...
##  $ magrange       : Factor w/ 3 levels "Weak[4-5)","Moderate[5-6)",..: 1 1 1 2 1 1 2 2 1 1 ...

Note: I attempted to color the earthquake types by magnitude range (weak, moderate, strong) by creating a color palette with color_factor but I failed.

Map Design

Choosing the correct icon for the leaflet interactive map is important to help the reader interpret the data correctly. Each of the red points below are single earthquakes and enhances the impact of the visualization with a red for alert color choice.

#* Create the map object and add circle marker*
leaflet(data = quakes1) %>% #*Initialize leaflet, pass on the quakes1 data*
  addProviderTiles("Esri.WorldImagery") %>% #* Add the third party tile*
  addCircleMarkers(lng = ~ longitude, lat = ~ latitude, #* Add circle earthquakes, long/lat*
                   radius = 4,
                   color = "red",
                   label = paste("Magnitude=", quakes1$mag, "Type=", quakes1$magrange))

Figure 1. World earthquakes greater than 4.5 in the past 30 days. Most earthquakes occurred along the archipelago coastal regions of Southeast Asia and along the west coast of South America with a focus on Chile. The earthquake icons have been enabled by hovering to show earthquake type (by magnitude) weak, moderate, strong as well as magnitude type.

I decided to zoom in to the west coast of Chile near Santiago as this area is typically one of the hardest hit.

  #### *Pass these colors to our earthquake point data (circles)*
#### *Label with magnitude strength when hovering over points*

leaflet() %>% 
   setView(lng=-66.1653, lat = -17.4139, zoom = 5) %>%
  addProviderTiles(providers$Esri.WorldImagery) %>%
  
    addCircleMarkers(data = quakes_greater4.5, label = ~place, radius = 3, color = "red") 
## Assuming "longitude" and "latitude" are longitude and latitude, respectively

Figure 2. Earthquakes greater than 4.5 in the past 30 days in southern Peru, northern Chile and Argentina. The earthquake icons have been enabled by hovering to show earthquakes labeled with place. Today June 12, 2023 Cochabamba had an earthquake 4.8 magnitude.

Cochabamba, Earthquake June 12, 2023 at 12:37pm magnitude 4.8.

Where in the world are the deepest quakes greater than 4.5 magnitude in the past 30 days?

#### Add a Legend to enable reader to recall what the colored values mean
#### Color the quakes by magnitude strength

quakesgreater6.0 <- quakes_greater4.5 %>%
  filter(mag >= 6.0) # create data frame with quakes than 6.0 magnitude

#### Call the colorNumeric() function to create a new palette function

pal <- colorNumeric(
  palette = "Reds",
  domain = quakesgreater6.0$depth)
#### Add a Legend to enable reader to recall what the colored values mean
#### Color the quakes by magnitude strength

map1 <- leaflet(data = quakes_greater4.5) %>% # call leaflet function and add data
  addProviderTiles(providers$Esri.WorldImagery) %>% # add tiles
  addCircleMarkers(color = ~pal(depth),
                   radius = 5,
                   fillOpacity = .5,
                   stroke = FALSE,
                   popup = paste0("mag"),
                   label = ~place) # add earthquake points
## Assuming "longitude" and "latitude" are longitude and latitude, respectively
## Warning in pal(depth): Some values were outside the color scale and will be
## treated as NA

## Warning in pal(depth): Some values were outside the color scale and will be
## treated as NA
title = "World Earthquakes greater than 4.5 magnitude in past 30 days"

map1

Figure 3. Earthquakes greater than 4.5 in the past 30 days colored by depth(km) and icons are labeled by place name. A deep earthquake occurred 2.0km southeast of Canilla, Guatemala.

Conclusion

Future Work

-Plot high risk earthquakes with high population areas and low income to assess which areas are most at risk for injury or worse death.

References

[USGS earthquake data past month] (https://earthquake.usgs.gov/earthquakes/feed/v1.0/csv.php)