Earthquake Magnitude / Depth Chart

[This article was first published on Exegetic Analytics » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I am working on a project related to secondary effects of earthquakes. To guide me in the analysis I need a chart showing the location, magnitude and depth of recent earthquakes. There are a host of such charts available already, but since I had the required data on hand, it seemed like a good idea to take a stab at it myself.

Getting the Data

The data was sourced from the US Geological Survey web site. I selected dates for the decade between between 1 January 2004 and 1 January 2013, magnitudes greater than 5 and chose CSV as the output format.

Loading the data into R is then simple. Some small transformations are required in order to interpret the time field in the data. I discarded a few columns which were not going to be useful, and added fields for the year and date of observation (for convenience alone: these data were already in the time field).

> catalog <- read.csv(file.path("data", "earthquake-catalog.csv"), stringsAsFactors = FALSE)
> #
> catalog <- within(catalog, {
+   time <- sub("T", " ", time)
+   time <- sub("Z", "", time)
+   time <- strptime(time, format = "%Y-%m-%d %H:%M:%S")
+   date <- as.Date(time)
+   year <- as.integer(strftime(time, format = "%Y"))
+ })
>
> catalog <- catalog[, c(12, 16, 17, 1, 2:5, 14)]

This is what the resulting data frame looks like:

> head(catalog)
          id year       date                time latitude longitude depth mag                           place
1 usc000lv53 2013 2013-12-31 2013-12-31 23:41:47  19.1673  120.0807 10.28 5.2  92km NW of Davila, Philippines
2 usc000lv0r 2013 2013-12-31 2013-12-31 21:32:01  19.1223  120.1797 10.00 5.2 83km NNW of Davila, Philippines
3 usb000m2uh 2013 2013-12-31 2013-12-31 20:04:32  19.0589  120.3057 20.61 5.0 70km NNW of Davila, Philippines
4 usc000luwe 2013 2013-12-31 2013-12-31 20:01:06  19.1181  120.2719 10.95 5.7 77km NNW of Burgos, Philippines
5 usb000m2ub 2013 2013-12-31 2013-12-31 13:55:02 -17.6528 -173.6869 15.38 5.0      114km NNE of Neiafu, Tonga
6 usc000lumu 2013 2013-12-31 2013-12-31 08:36:30 -15.6555 -172.9340 31.82 5.1       93km ENE of Hihifo, Tonga

Making the Charts

Time to generate those charts. There are lots of ways to make maps in R, I chose to use a generic option: ggplot2.

> require(ggplot2)
> require(maps)
> require(grid)
> 
> world.map <- map_data("world")
> 
> ggplot() +
+   geom_polygon(data = world.map, aes(x = long, y = lat, group = group),
+                fill = "#EEEECC") +
+   geom_point(data = catalog, alpha = 0.25,
+              aes(x = longitude, y = latitude, size = mag, colour = depth)) +
+   labs(x = NULL, y = NULL) +
+   scale_colour_gradient("Depth [m]", high = "red") +
+   scale_size("Magnitude") +
+   coord_fixed(ylim = c(-82.5, 87.5), xlim = c(-185, 185)) +
+   theme_classic() +
+   theme(axis.line = element_blank(), axis.text = element_blank(),
+         axis.ticks = element_blank(),
+         plot.margin=unit(c(3, 0, 0, 0),"mm"),
+         legend.text = element_text(size = 6),
+         legend.title = element_text(size = 8, face = "plain"),
+         panel.background = element_rect(fill='#D6E7EF'))

The resulting plot gives the location of the earthquakes as points, with magnitudes indicated by the sizes of the points and depths given by their colour.

earthquake-map

The Earth’s tectonic plates are well defined by the numerous interplate earthquakes, and there is a liberal sprinkling of intreplate events as well.

I made another chart showing the distribution of earthquakes broken down by year.

earthquake-map-panels

Distribution of Earthquake Magnitudes

While we are taking a high level look at the data, it’s interesting to see how the magnitudes are distributed. A logartihmic scale is necessary to make the frequencies visible over the full range of magnitudes.

ggplot(catalog, aes(x = mag)) +
  xlab("Magnitude") + ylab("Number of Earthquakes") +
  stat_bin(drop = TRUE, binwidth = 0.25) +
  scale_y_log10(breaks = c(1, 10, 100, 1000)) +
  theme_classic()

earthquake-magnitude-histogram

Very nice: consistent with a power law, as described by the Gutenberg–Richter law.

To leave a comment for the author, please follow the link and comment on their blog: Exegetic Analytics » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)