1. Data.

Get some earthquake data: all M6.0 or greater quakes since 1980 as recorded by the USGS.

world_1980_m6 <- read.csv('world_1980_m6.csv')
world_1980_m6$time2 <- strftime(world_1980_m6$time,
                  format = "%Y")
start.time <- strftime('1980-01-01', format="%Y-%m-%d")
world_1980_m6$ndays <- as.numeric(difftime(strftime(world_1980_m6$time, format="%Y-%m-%d"),
                                           start.time, units = "days"))
head(world_1980_m6)

2. Magnitude and Depth

Plot a histogram of magnitudes.

library(ggplot2)
require(gridExtra)
plot1 <- qplot(world_1980_m6$mag, geom="histogram",
               main = "Magnitude Distribution",
               xlab = "Magnitude")
plot2 <- qplot(world_1980_m6$depth, geom="histogram",
               main = "Depth Distribution",
               xlab = "Depth")
grid.arrange(plot1, plot2, ncol=2)

Magnitude and depth both exhibit long right tails. Let’s examine the relationship further by looking at k-means clustering of magnitude-depth pairs.

2.1Mag-Depth K-means

Most earthquakes, even above M6.0, occur at a depth of 200 or less. There is a small break where the earth appears to be more stable around 300 depth. However, earthquakes begin occuring again with regularity between 500 and 600 depth.

set.seed(20)
world_1980_m6.cluster <- kmeans(world_1980_m6[, 4:5], 2, nstart = 20)
world_1980_m6$cluster <- as.factor(world_1980_m6.cluster$cluster)
ggplot(world_1980_m6, aes(depth, mag, color = world_1980_m6$cluster)) + geom_point()

2.2 Mag-Depth Kernel Density

Are deeper earthquakes stronger? We can use our clusters to visualize magnitude distribution for each depth range.

world_1980_m6 %>%
  ggvis(~mag, fill = ~cluster) %>%
  group_by(cluster) %>%
  layer_densities()

3. Frequency

Plot as a time series.

library(dplyr)
library(ggvis)
world_1980_m6 %>% 
  group_by(time2) %>% 
  summarise(count = n()) %>% 
  ggvis(~time2, ~count) %>% 
  layer_bars() %>% 
  add_axis("x",
           values = seq(1980, 2017, by = 5),
           properties = axis_props(
             labels = list(angle = 90, align = "left", fontSize = 20)
           ),
           tick_size_major = 10,
           tick_size_minor = 5
  )

3.1 Frequency and Magnitude

world_1980_m6 %>% 
  ggvis(~ndays, ~mag) %>% 
  layer_points()
world_1980_m6 %>% 
  ggvis(~time2, ~mag) %>% 
  layer_boxplots() %>% 
  add_axis("x",
           values = seq(1980, 2017, by = 5),
           properties = axis_props(
             labels = list(angle = 90, align = "left", fontSize = 20)
           ),
           tick_size_major = 10,
           tick_size_minor = 5
  )

NA

4. Geographic Distribution

library(leaflet)
pal <- colorNumeric("Reds", c(6, 9), na.color = "#808080")
leaflet(world_1980_m6) %>% 
  addTiles() %>% 
  addCircleMarkers(~longitude,
                   ~latitude,
                   popup = ~place,
                   fillOpacity = 0.9,
                   clusterOptions = markerClusterOptions())
LS0tCnRpdGxlOiAiU29tZSBFYXJ0aHF1YWtlcyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMgMS4gRGF0YS4KR2V0IHNvbWUgZWFydGhxdWFrZSBkYXRhOiBhbGwgTTYuMCBvciBncmVhdGVyIHF1YWtlcyBzaW5jZSAxOTgwIGFzIHJlY29yZGVkIGJ5IHRoZSBVU0dTLgoKYGBge3J9CndvcmxkXzE5ODBfbTYgPC0gcmVhZC5jc3YoJ3dvcmxkXzE5ODBfbTYuY3N2JykKd29ybGRfMTk4MF9tNiR0aW1lMiA8LSBzdHJmdGltZSh3b3JsZF8xOTgwX202JHRpbWUsCiAgICAgICAgICAgICAgICAgIGZvcm1hdCA9ICIlWSIpCnN0YXJ0LnRpbWUgPC0gc3RyZnRpbWUoJzE5ODAtMDEtMDEnLCBmb3JtYXQ9IiVZLSVtLSVkIikKd29ybGRfMTk4MF9tNiRuZGF5cyA8LSBhcy5udW1lcmljKGRpZmZ0aW1lKHN0cmZ0aW1lKHdvcmxkXzE5ODBfbTYkdGltZSwgZm9ybWF0PSIlWS0lbS0lZCIpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RhcnQudGltZSwgdW5pdHMgPSAiZGF5cyIpKQpoZWFkKHdvcmxkXzE5ODBfbTYpCmBgYAoKIyMgMi4gTWFnbml0dWRlIGFuZCBEZXB0aApQbG90IGEgaGlzdG9ncmFtIG9mIG1hZ25pdHVkZXMuCgpgYGB7ciBlcS1tYWctaGlzdH0KbGlicmFyeShnZ3Bsb3QyKQpyZXF1aXJlKGdyaWRFeHRyYSkKcGxvdDEgPC0gcXBsb3Qod29ybGRfMTk4MF9tNiRtYWcsIGdlb209Imhpc3RvZ3JhbSIsCiAgICAgICAgICAgICAgIG1haW4gPSAiTWFnbml0dWRlIERpc3RyaWJ1dGlvbiIsCiAgICAgICAgICAgICAgIHhsYWIgPSAiTWFnbml0dWRlIikKcGxvdDIgPC0gcXBsb3Qod29ybGRfMTk4MF9tNiRkZXB0aCwgZ2VvbT0iaGlzdG9ncmFtIiwKICAgICAgICAgICAgICAgbWFpbiA9ICJEZXB0aCBEaXN0cmlidXRpb24iLAogICAgICAgICAgICAgICB4bGFiID0gIkRlcHRoIikKZ3JpZC5hcnJhbmdlKHBsb3QxLCBwbG90MiwgbmNvbD0yKQpgYGAKCk1hZ25pdHVkZSBhbmQgZGVwdGggYm90aCBleGhpYml0IGxvbmcgcmlnaHQgdGFpbHMuIExldCdzIGV4YW1pbmUgdGhlIHJlbGF0aW9uc2hpcCBmdXJ0aGVyIGJ5IGxvb2tpbmcgYXQgay1tZWFucyBjbHVzdGVyaW5nIG9mIG1hZ25pdHVkZS1kZXB0aCBwYWlycy4KCiMjIyMgMi4xTWFnLURlcHRoIEstbWVhbnMKCk1vc3QgZWFydGhxdWFrZXMsIGV2ZW4gYWJvdmUgTTYuMCwgb2NjdXIgYXQgYSBkZXB0aCBvZiAyMDAgb3IgbGVzcy4gVGhlcmUgaXMgYSBzbWFsbCBicmVhayB3aGVyZSB0aGUgZWFydGggYXBwZWFycyB0byBiZSBtb3JlIHN0YWJsZSBhcm91bmQgMzAwIGRlcHRoLiBIb3dldmVyLCBlYXJ0aHF1YWtlcyBiZWdpbiBvY2N1cmluZyBhZ2FpbiB3aXRoIHJlZ3VsYXJpdHkgYmV0d2VlbiA1MDAgYW5kIDYwMCBkZXB0aC4KYGBge3J9CnNldC5zZWVkKDIwKQoKd29ybGRfMTk4MF9tNi5jbHVzdGVyIDwtIGttZWFucyh3b3JsZF8xOTgwX202WywgNDo1XSwgMiwgbnN0YXJ0ID0gMjApCgp3b3JsZF8xOTgwX202JGNsdXN0ZXIgPC0gYXMuZmFjdG9yKHdvcmxkXzE5ODBfbTYuY2x1c3RlciRjbHVzdGVyKQoKZ2dwbG90KHdvcmxkXzE5ODBfbTYsIGFlcyhkZXB0aCwgbWFnLCBjb2xvciA9IHdvcmxkXzE5ODBfbTYkY2x1c3RlcikpICsgZ2VvbV9wb2ludCgpCgpgYGAKCiMjIyMgMi4yIE1hZy1EZXB0aCBLZXJuZWwgRGVuc2l0eQpBcmUgZGVlcGVyIGVhcnRocXVha2VzIHN0cm9uZ2VyPyBXZSBjYW4gdXNlIG91ciBjbHVzdGVycyB0byB2aXN1YWxpemUgbWFnbml0dWRlIGRpc3RyaWJ1dGlvbiBmb3IgZWFjaCBkZXB0aCByYW5nZS4KCmBgYHtyfQp3b3JsZF8xOTgwX202ICU+JQogIGdndmlzKH5tYWcsIGZpbGwgPSB+Y2x1c3RlcikgJT4lCiAgZ3JvdXBfYnkoY2x1c3RlcikgJT4lCiAgbGF5ZXJfZGVuc2l0aWVzKCkKYGBgCgojIyAzLiBGcmVxdWVuY3kKUGxvdCBhcyBhIHRpbWUgc2VyaWVzLgoKYGBge3IgZXEtdGltZXNlcmllc30KbGlicmFyeShkcGx5cikKbGlicmFyeShnZ3ZpcykKCndvcmxkXzE5ODBfbTYgJT4lIAogIGdyb3VwX2J5KHRpbWUyKSAlPiUgCiAgc3VtbWFyaXNlKGNvdW50ID0gbigpKSAlPiUgCiAgZ2d2aXMofnRpbWUyLCB+Y291bnQpICU+JSAKICBsYXllcl9iYXJzKCkgJT4lIAogIGFkZF9heGlzKCJ4IiwKICAgICAgICAgICB2YWx1ZXMgPSBzZXEoMTk4MCwgMjAxNywgYnkgPSA1KSwKICAgICAgICAgICBwcm9wZXJ0aWVzID0gYXhpc19wcm9wcygKICAgICAgICAgICAgIGxhYmVscyA9IGxpc3QoYW5nbGUgPSA5MCwgYWxpZ24gPSAibGVmdCIsIGZvbnRTaXplID0gMjApCiAgICAgICAgICAgKSwKICAgICAgICAgICB0aWNrX3NpemVfbWFqb3IgPSAxMCwKICAgICAgICAgICB0aWNrX3NpemVfbWlub3IgPSA1CiAgKQpgYGAKCiMjIyMgMy4xIEZyZXF1ZW5jeSBhbmQgTWFnbml0dWRlCmBgYHtyIGZyZXEtbWFnLWRvdH0Kd29ybGRfMTk4MF9tNiAlPiUgCiAgZ2d2aXMofm5kYXlzLCB+bWFnKSAlPiUgCiAgbGF5ZXJfcG9pbnRzKCkKYGBgCgpgYGB7ciBmcmVxLW1hZy1ib3h9CndvcmxkXzE5ODBfbTYgJT4lIAogIGdndmlzKH50aW1lMiwgfm1hZykgJT4lIAogIGxheWVyX2JveHBsb3RzKCkgJT4lIAogIGFkZF9heGlzKCJ4IiwKICAgICAgICAgICB2YWx1ZXMgPSBzZXEoMTk4MCwgMjAxNywgYnkgPSA1KSwKICAgICAgICAgICBwcm9wZXJ0aWVzID0gYXhpc19wcm9wcygKICAgICAgICAgICAgIGxhYmVscyA9IGxpc3QoYW5nbGUgPSA5MCwgYWxpZ24gPSAibGVmdCIsIGZvbnRTaXplID0gMjApCiAgICAgICAgICAgKSwKICAgICAgICAgICB0aWNrX3NpemVfbWFqb3IgPSAxMCwKICAgICAgICAgICB0aWNrX3NpemVfbWlub3IgPSA1CiAgKQogIApgYGAKCiMjIDQuIEdlb2dyYXBoaWMgRGlzdHJpYnV0aW9uCmBgYHtyIGVxLW1hcH0KbGlicmFyeShsZWFmbGV0KQpwYWwgPC0gY29sb3JOdW1lcmljKCJSZWRzIiwgYyg2LCA5KSwgbmEuY29sb3IgPSAiIzgwODA4MCIpCmxlYWZsZXQod29ybGRfMTk4MF9tNikgJT4lIAogIGFkZFRpbGVzKCkgJT4lIAogIGFkZENpcmNsZU1hcmtlcnMofmxvbmdpdHVkZSwKICAgICAgICAgICAgICAgICAgIH5sYXRpdHVkZSwKICAgICAgICAgICAgICAgICAgIHBvcHVwID0gfnBsYWNlLAogICAgICAgICAgICAgICAgICAgZmlsbE9wYWNpdHkgPSAwLjksCiAgICAgICAgICAgICAgICAgICBjbHVzdGVyT3B0aW9ucyA9IG1hcmtlckNsdXN0ZXJPcHRpb25zKCkpCmBgYA==