library(leaflet)
library(sf)
library(tidyverse)
library(knitr)Japan Earthquakes Geo Tutorial
Japan Earthquakes
Geo mapping using Japan Earthquakes 2001-2018 dataset. Source: USGS Earthquake Catalog (https://earthquake.usgs.gov/earthquakes/search/).
This tutorial is adapted from: https://towardsdatascience.com/how-to-make-stunning-geomaps-in-r-a-complete-guide-with-leaflet-be1b857f1644
This is a large dataset with over 14,000 observations and 22 variables, including:
| Variables Used |
|---|
| time |
| longitude |
| latitude |
| depth |
| mag (magnitude) |
| magType (magnitude type) |
| place |
| type |
| horizontalError |
| depthError |
| magError |
Load the necessary libraries
Set the lat and long values for Japan
japan_lon <- 138.129731
japan_lat <- 36.2615855Set the working directory and read in the earthquake data
quakes <- read_csv("japan_quakes_01-18.csv")Look at Japan’s earthquakes over time
First create a new variable date that removes the time stamp from the time variable
library(lubridate)
quakes$date <- as.Date(format(ymd_hms(quakes$time),
format = '%Y-%m-%d'))
head(quakes)# A tibble: 6 × 23
time latitude longitude depth mag magType nst gap dmin
<dttm> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
1 2018-11-27 14:34:20 48.4 155. 35 4.9 mb NA 92 5.04
2 2018-11-26 23:33:50 36.1 140. 48.8 4.8 mww NA 113 1.36
3 2018-11-26 13:04:02 38.9 142. 50.6 4.5 mb NA 145 1.29
4 2018-11-26 05:20:16 50.1 156. 66.3 4.6 mb NA 128 3.19
5 2018-11-25 09:19:05 34.0 134. 38.2 4.6 mb NA 104 0.558
6 2018-11-25 03:16:46 48.4 155. 35 4.6 mb NA 85 4.99
# ℹ 14 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>, date <date>unique(quakes$magType)[1] "mb" "mww" "mwr" "mwb" "mwc" "ms" "mw" "m" Create a scatterplot graph of earthquakes with magnitudes greater than 5
USGS classifies different types of magnitudes. You can find more information here: https://www.usgs.gov/programs/earthquake-hazards/magnitude-types
quakes1 <- quakes |>
filter(mag >=5) |>
group_by(date) |>
reframe(avg_mag = mean(mag), magType)
ggplot(quakes1, aes(x=date, y=avg_mag, color = magType)) +
geom_point(alpha = 0.4) +
labs(title = "Earthquakes in Japan by Magnitude Type",
caption = "Source: USGS") +
facet_wrap(~magType) +
scale_color_viridis_d("Magnitude Type") +
theme_bw()
Here, we can see that mww and mb types were more prevalent after 2011 and prior to 2011, mwc and mwb were more prevalent types.
Create a scatterplot of magnitude versus depth
ggplot(quakes,
aes(x=depth,
y=mag,
color = magType)) +
geom_point(alpha = 0.1) +
scale_color_viridis_d()+
geom_jitter() +
labs(title = "Earthquakes in Japan by Magnitude Type",
caption = "Source: USGS") +
theme_bw()
This plot appears to have two clumps with peaks at <50 meters depth and about 400-500 meters depth
Try to Facet for the magType
ggplot(quakes, aes(x=depth,
y=mag,
color = magType)) +
geom_point(alpha = 0.05) +
scale_color_viridis_d()+
geom_jitter() +
facet_wrap(~magType) +
labs(title = "Earthquakes in Japan by Magnitude Type",
caption = "Source: USGS") +
theme_bw()
Since there are so few, remove magTypes: m, ms, and mw
quakes2 <- quakes |>
filter(magType %in% c("mb",
"mwb",
"mwc",
"mwr",
"mww"))
ggplot(quakes2, aes(x=depth,
y=mag,
color = magType)) +
geom_point(alpha = 0.05) +
scale_color_viridis_d()+
geom_jitter() +
facet_wrap(~magType) +
labs(title = "Earthquakes in Japan by Magnitude Type",
caption = "Source: USGS") +
theme_bw()
Mwc, mwb, and mww in particular seem to increase in magnitude for depths greater than 400 meters.
Explore the duration
ggplot(quakes,
aes(dmin)) +
geom_density(bins = 15)Warning in geom_density(bins = 15): Ignoring unknown parameters: `bins`Warning: Removed 10485 rows containing non-finite outside the scale range
(`stat_density()`).
Most earthquakes last for under 7 minutes
Explore the outliers
ggplot(quakes2,
aes(x=dmin,
color = magType))+
geom_density(alpha = 0.4) +
scale_color_viridis_d() +
theme_bw() +
labs(title = "Earthquakes in Japan by Depth",
caption = "Source: USGS") Warning: Removed 10444 rows containing non-finite outside the scale range
(`stat_density()`).
mww has some of the longest durations with the greatest magnitudes and depths
quakes_mww <- quakes |>
filter(magType == "mww") |>
filter(!is.na(mag) &
!is.na(depth) &
!is.na(dmin))
ggplot(quakes_mww,
aes(x=date,
y=mag,
color = dmin)) +
geom_point(aes(size=depth)) +
scale_color_gradient(high = "#14010d",
low = "#f2079c") +
theme_bw() +
labs(title = "MWW Earthquakes in Japan by Over the Years Based \n on Magnitude, Depth, and Duration",
caption = "Source: USGS") 
Filter for magnitudes greater than 6
strong <- quakes |>
filter(mag >= 6)Mapping Earthquakes
Decide what style you would like your map
Use this link to look at all the options: https://leaflet-extras.github.io/leaflet-providers/preview/
Draw a first map using leaflet and the Esri World Street Map
Use the function addProviderTiles() to decide which style you want - I have included a few different types for you to try out.
Some examples:
addProviderTiles(“Stamen.Terrain”)
addProviderTiles(“Stamen.Watercolor”)
addProviderTiles(“Esri.WorldPhysical”)
addProviderTiles(“Esri.NatGeoWorldMap”)
First, we need to set the location of the map of Japan. A quick google search gives 36.2 degrees North, 138.2 degrees East. That translates to lng = 138.2, lat = 36.2
(North = + lat, South = - lat, East = + lng, West = - lng)
Zoom provides the level of granularity. Play with zoom = 1 versus zoom = 6 to see what the initial setting is.
leaflet() |>
setView(lng = 138.2,
lat = 36.2,
zoom =6) |>
addProviderTiles("Esri.WorldStreetMap") |>
addCircles(data = strong,
radius = strong$mag)Assuming "longitude" and "latitude" are longitude and latitude, respectivelyTweak the marker size
You can use the following completely non-scientific formula to calculate marker size:
The x here represents the magnitude, and c represents a constant you can play around with. The larger it is, the easier it is to spot stronger earthquakes.
Implement this formula in the radius with the value for c is set to 2.
leaflet() |>
setView(lng = japan_lon,
lat = japan_lat,
zoom = 6) |>
addProviderTiles("Esri.WorldStreetMap") |>
addCircles(data = strong,
radius = sqrt(10^strong$mag)*2)Assuming "longitude" and "latitude" are longitude and latitude, respectivelyTweak the markers
You can tweak the markers using addCircles() in terms of color, fillColor, and fillOpacity
leaflet() |>
setView(lng = japan_lon,
lat = japan_lat,
zoom = 6) |>
addProviderTiles("Esri.WorldStreetMap") |>
addCircles(data = strong,
radius = sqrt(10^strong$mag) * 2,
color = "#14010d",
fillColor = "#f2079c",
fillOpacity = 0.25)Assuming "longitude" and "latitude" are longitude and latitude, respectivelyCreate a popup using paste0
create a line break using < br >
surround text with < b > makes it bold
popupquake <- paste0(
"<b>Time: </b>", strong$time, "<br>",
"<b>Magnitude: </b>", strong$mag, "<br>",
"<b>Depth (km): </b>", strong$depth, "<br>",
"<strong>Place: </strong>", strong$place, "<br>")Add the popup to the map
Click on the points to see the popup tooltip for details about each point
leaflet() |>
setView(lng = japan_lon,
lat = japan_lat,
zoom = 6) |>
addProviderTiles("Esri.WorldStreetMap") |>
addCircles(
data = strong,
radius = sqrt(10^strong$mag) * 2,
color = "#14010d",
fillColor = "#f2079c",
fillOpacity = 0.35,
popup = popupquake)Notice the largest point is in 2011, with a magnitude of 9.1, which caused the Fukushima nuclear power plant meltdown.