A Byte of my 2.2-lb Brain

Just sharing stuff…

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.bohol earthquake.png

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")

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Information

This entry was posted on January 18, 2017 by in Blog, Geek, Philippines and tagged , , , , , , .
%d bloggers like this: