Earthquake map (Greece)

Earthquakes: Magnitude / Depth Chart

The following code produces a chart showing the location, magnitude and depth of recent earthquakes in Greece. 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.

The data was sourced from the US Geological Survey web site. (http://earthquake.usgs.gov/earthquakes/search/).

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.

suppressPackageStartupMessages(library(rgdal))
suppressPackageStartupMessages(library(maptools))
suppressPackageStartupMessages(require(grid))
#library(gpclib)
suppressPackageStartupMessages(library(ggthemes))

setwd("D:/GD/R_data/statistics")

catalog <- read.csv("../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)]


catalog$year <-  floor(catalog$year/10)*10


#download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries.zip", "../data/countries.zip")
# Then unzip
#unzip("../data/countries.zip",exdir="../data")


#gpclibPermit()
world.map <- readOGR(dsn="../data", layer="ne_10m_admin_0_countries")

## OGR data source with driver: ESRI Shapefile 
## Source: "../data", layer: "ne_10m_admin_0_countries"
## with 255 features and 65 fields
## Feature type: wkbPolygon with 2 dimensions

world.ggmap <- fortify(world.map, region = "NAME")

## Loading required package: rgeos
## rgeos version: 0.3-6, (SVN revision 450)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE

world.ggmap <- world.ggmap[world.ggmap$id=="Greece",]

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.

Distribution of Earthquake Magnitudes

It’s interesting to see how the magnitudes are distributed. The greener the color, the larger the focal depth.

ggplot(catalog, aes(x = mag, fill=factor(depth))) +  
  scale_fill_discrete(h = c(0, 130) + 15, c = 100, l = 65, guide="none") +
  xlab("Magnitude") + ylab("Number of Earthquakes") +
  geom_histogram(binwidth = 0.25) +
  theme_stata() +
  ggtitle("Distribution of Earthquake Magnitudes")+ 
  theme(legend.position = "bottom", 
        axis.text.y=element_text(angle=0,hjust=1.3))