Earthquake Maps in R

I am publishing the R code I wrote to create the map of Bohol shown below with earthquake location points.

In October 2013, a very strong earthquake, with a 7.2 magnitude, hit the island province Bohol, which is located in Central Visayas, Philippines. According to reports, the natural disaster claimed at least 222 lives on 15 October 2013; around 950 were injured. It’s the strongest and deadliest earthquake that struck the archipelago in recent years.  Just three weeks after the quake, Super Typhoon Haiyan hit the very same Central Visayas region. 😦

On 18 October 2013, I decided to map the earthquake locations within the Bohol region from 15-17 October 2013. The red data point on the map is the location of the magnitude 7.2 earthquake that hit Bohol on 15 October 2013.

Providing the R code below as is and without comments. If you want a step-by-step guide to the process, check the R markdown document I wrote here; it’s pretty verbose. The document also includes  info on the data sources (for both the map and USGS earthquake data).

library(maptools)

gadm<-readRDS("../PH_GIS_Data/PHL_adm2.rds")
bohol<-gadm[gadm$NAME_1=="Bohol",]
catigbian<-bohol[bohol$NAME_2=="Catigbian",]

quake<-read.csv('equake3.csv')
quake<-quake[quake$country=="Philippines",]
quake$date <- as.Date(quake$date, format="%d/%m/%y")
quake <- quake[quake$date >= as.Date("2013-10-15") & quake$date <= as.Date("2013-10-17"),]

malakas = quake[quake$mag==max(quake$mag),]

colors <- c("salmon", "deepskyblue2", "darkseagreen")

plot(bohol, axes=T, border='gray50', bg='gray80', col='black', lwd=1.3, xlim=c(123.68114,124.64547), ylim=c(9.48583, 10.4))
par(new=TRUE)
plot(catigbian, col="yellow", xlim=c(123.68114,124.64547), ylim=c(9.48583, 10.4), bg="transparent")

points(quake$longitude,quake$latitude, col=colors[quake$date], pch=16, cex = sqrt(quake$mag/pi))
points(malakas$longitude,malakas$latitude, col='red', pch=16, cex = sqrt(malakas$mag/2))
points(malakas$longitude,malakas$latitude, col='red', pch=1, cex = sqrt(malakas$mag))
points(malakas$longitude,malakas$latitude, col='red', pch=1, cex = 1.2*sqrt(malakas$mag))
legend(124.5,10.36, legend = levels(quake$date), col=colors, cex = 0.9, pch = 16, bty = "n")