setwd("/home/sam/Documents/R")
library(mapdata)

#Declare neg() function
neg<-function(x) -x

#Read in data
GNS<-read.csv("NZ_earthquakes2010.csv") #Earthquake data
loc<-read.csv("Canterbury_towns.csv") #Canterbury gazetter
loc<-loc[-6,]

#Correct formatting
#Date
GNS$Date<-as.POSIXct(paste(GNS$ORI_YEAR, GNS$ORI_MONTH, GNS$ORI_DAY, GNS$ORI_HOUR, GNS$ORI_MINUTE, GNS$ORI_SECOND, sep=" "), format="%Y %m %d %H %M %S", tz="GMT")

GNS$LONG_cor<-if(GNS$LONG< 0) GNS$LONG+180 else GNS$LONG-180


#Gazetter lat and longs
loc$LAT<-loc$Lat_deg-loc$Lat_min/60
loc$LONG<-loc$Long_deg+loc$Long_min/60


#Subset all earthquakes in Canterbury
cant<-GNS[GNS$LAT<neg(42.4) & GNS$LAT>neg(44) & GNS$LONG> 171 & GNS$LONG< 173.5,]
cant2<-cant[-1,] #One outlier distance-wise


#Subset all Canterbury earthquakes before and after the big one.
after<-cant[cant$Date>as.POSIXct("2010-09-03"),]
before<-cant[cant$Date<as.POSIXct("2010-09-03"),]

#New variables
after$time_from<-after$Date-after$Date[2]   #Time afterward
after$time_diff<-rep(NA, length(after$Date))
for(i in 2:length(after$Date)) after$time_diff[i]<-after$Date[i]-after$Date[i-1]

after$dist<-as.matrix(dist(cbind(cant$LAT,cant$LONG)))[,2] #Distance from epicentre of the first

#Number of aftershocks
dim(after)

#Number of aftershocks of greater than 4, 4.5 and 5 magnitude
mags<-c(4,4.5,5)
sapply(mags, FUN=function(x)length(which(after$MAG>x)))

bigger<-lapply(mags, FUN=function(x)after[after$MAG>x,])

biggerb<-lapply(mags, FUN=function(x)before[before$MAG>x,])

vec<-NULL
for(i in 1:(length(cant$MAG)-4)) vec[i]<-mean(cant$MAG[i:(i+4)])

#######################
#Plots

#Location of 2010 NZ earthquakes
plot(GNS$LAT~GNS$LONG_cor, type="n")
points(GNS$LAT~GNS$LONG_cor, cex=GNS$MAG)

map("nzHires", xlim=c(165,180), ylim=c(-48,-34))
map.axes()
mtext("2010 Earthquakes of magnitude less than 4\n")
points(GNS$LAT[GNS$MAG<4]~GNS$LONG[GNS$MAG<4], pch=19)

map("nzHires", xlim=c(165,180), ylim=c(-48,-34))
map.axes()
mtext("2010 Earthquakes of magnitude greater than 4\n")
points(GNS$LAT[GNS$MAG>4]~GNS$LONG[GNS$MAG>4], pch=19, col="orange")
points(GNS$LAT[GNS$MAG>5]~GNS$LONG[GNS$MAG>5], pch=19, col="blue")
legend(166, -34, legend=c("> 4", "> 5"), pch=19, col=c("orange", "blue"), title="Magnitude")

#Distribution of Magnitude and Depth
hist(GNS$MAG, main="Magnitude of New Zealand earthqukes \n detected in 2010 to 8 September", xlab="Magnitude (Richter scale)", col="grey25")

hist(GNS$DEPTH, main="Depth of New Zealand earthquakes \ndetected in 2010 to 8 September", xlab="Depth (km)", col="grey25")


#Magnitude related to date
#All New Zealand
plot(GNS$MAG ~ GNS$Date, type="p", pch=18, ylab="Magnitude", xlab="Date")
points(cant$MAG ~ cant$Date, pch=18, col="blue")
legend(1265599232, 7, legend=c("Canterbury", "Elsewherein NZ"), pch=18, col=c("blue", "black"))

#Just aftershocks, with 
plot(after$MAG~after$Date, type="l", xlab="Time", ylab="Magnitude")
lines(vec ~ cant$Date[1:length(vec)], col="blue", lwd=2)

plot(GNS$MAG ~GNS$DEPTH, pch=19, xlab="Depth (km)", ylab="Magnitude")
points(after$MAG ~ after$DEPTH, pch=19, col="blue")


plot(after$MAG ~ after$time_from)

plot(after$DEPTH ~ after$time_from)

plot(after$time_diff ~ after$Date)

plot(after$MAG ~ after$dist)

plot(after$DEPTH ~ after$dist)

plot(after$MAG ~ after$time_diff)


#-----------------------------------------
#Canterbury aftershocks prior to the earthquake

map("nzHires", regions="South Island", xlim=c(171,173.5), ylim=c(-44,-43.1), fill=TRUE, col="light green")
mtext("Earthquakes in the Canterbury region in 2010 prior to 4 Sep 2010 \n")
map.axes()
points(before$LAT ~ before$LONG, pch=19)
text(before$LAT[86] ~ before$LONG[86], label=paste(before$Date[86], "Mag:", before$MAG[86]), adj=c(0.5,-1))
text(loc$LONG, loc$LAT, labels=loc$Locality, col="blue", adj=c(-0.1,0.5))
points(loc$LAT ~ loc$LONG, pch=19, col="blue", lwd=2)

#-----------------------------------------
#Only the larger ones

map("nzHires", regions="South Island", xlim=c(171,173.5), ylim=c(-44,-43.1), fill=TRUE, col="light green")
mtext("Location and magnitude of aftershocks from the Canterbury Earthquake \n 4--7 Sep 2010 \n")
map.axes()
points(after$LAT ~ after$LONG, pch=19)
points(bigger[[1]]$LAT ~ bigger[[1]]$LONG, pch=21, cex=2, bg="red")
points(bigger[[2]]$LAT ~ bigger[[2]]$LONG, pch=21, cex=2, bg="yellow")
points(bigger[[3]]$LAT ~ bigger[[3]]$LONG, pch=21, cex=2, bg="pink")
points(cant$LAT[cant$MAG> 7] ~ cant$LONG[cant$MAG> 7], pch="X", cex=2)

text(loc$LONG, loc$LAT, labels=loc$Locality, col="grey25", adj=c(-0.1,0.5))
points(loc$LAT ~ loc$LONG, pch=19, col="grey25", lwd=2)

legend(x=171.1, y=-43.2, legend=c("< 4", "> 4", "> 4.5", "> 5", "Epicentre"), pch=c(19, 21, 21, 21, 88), pt.cex=c(1,rep(2,3),1.5), pt.bg=c(NA, "red", "yellow", "pink", "black", "black"), title="Magnitude")

#-----------------------------------------

points(cant$LAT ~ cant$LONG, pch=19)




#######################
#Printing out the plots

png(file="all.png", width=600, height=600)
plot(GNS$MAG ~ GNS$Date, type="p", pch=18, ylab="Magnitude", xlab="Date", main="Earthquakes originating in New Zealand from 31 Aug--6 Sep 2010")
points(cant$MAG ~ cant$Date, pch=18, col="blue")
legend(1283609257, 7, legend=c("Canterbury", "Elsewherein NZ"), pch=18, col=c("blue", "black"))
dev.off()

#----------------------------------

png(file="all2010.png", width=900, height=600)
plot(GNS$MAG ~ GNS$Date, type="p", pch=18, ylab="Magnitude", xlab="Date", main="New Zealand earthquakes in 2010 to 8 September")
points(cant$MAG ~ cant$Date, pch=18, col="blue")
legend(1265599232, 7, legend=c("Canterbury", "Elsewherein NZ"), pch=18, col=c("blue", "black"))
dev.off()

#----------------------------------

png(file="line.png", width=600, height=600)
plot(cant$MAG~cant$Date, type="l", xlab="Time", ylab="Magnitude", main="Canterbury earthquakes, 4--6 Sep 2010")
dev.off()

#----------------------------------

png(file="depth.png", width=600, height=600)
plot(GNS$MAG ~GNS$DEPTH, pch=19, xlab="Depth (km)", ylab="Magnitude", main="New Zealand Earthquakes, 31 Aug-6 Sep 2010")
points(cant$MAG ~ cant$DEPTH, pch=19, col="blue")
legend(200, 7, legend=c("Canterbury", "Elsewherein NZ"), pch=18, col=c("blue", "black"))
dev.off()

#----------------------------------

png(file="after.png", width=900, height=600)
map("nzHires", regions="South Island", xlim=c(171,173.5), ylim=c(-44,-43.1), fill=TRUE, col="light green")
mtext("Location and magnitude of aftershocks from the Canterbury Earthquake \n 4--7 Sep 2010 \n", cex=2)
map.axes()
points(after$LAT ~ after$LONG, pch=19)
points(bigger[[1]]$LAT ~ bigger[[1]]$LONG, pch=21, cex=2, bg="red")
points(bigger[[2]]$LAT ~ bigger[[2]]$LONG, pch=21, cex=2, bg="yellow")
points(bigger[[3]]$LAT ~ bigger[[3]]$LONG, pch=21, cex=2, bg="pink")
points(cant$LAT[cant$MAG> 7] ~ cant$LONG[cant$MAG> 7], pch="X", cex=2)

text(loc$LONG, loc$LAT, labels=loc$Locality, col="dark blue", adj=c(-0.1,0.5), cex=2)
points(loc$LAT ~ loc$LONG, pch=19, col="dark blue", lwd=2, cex=2)

legend(x=171.1, y=-43.2, legend=c("< 4", "> 4", "> 4.5", "> 5", "Epicentre"), pch=c(19, 21, 21, 21, 88), pt.cex=c(1,rep(2,3),1.5), pt.bg=c(NA, "red", "yellow", "pink", "black", "black"), title="Magnitude")
dev.off()

#----------------------------------

png(file="before.png", width=900, height=600)
map("nzHires", regions="South Island", xlim=c(171,173.5), ylim=c(-44,-43.1), fill=TRUE, col="light green")
mtext("Earthquakes in the Canterbury region in 2010 prior to 4 Sep 2010 \n", cex=2)
map.axes()
points(before$LAT ~ before$LONG, pch=19)
text(before$LAT[86] ~ before$LONG[86], label=paste(before$Date[86], "Mag:", before$MAG[86]), adj=c(0.5,-1))
text(loc$LONG, loc$LAT, labels=loc$Locality, col="blue", adj=c(-0.1,0.5))
points(loc$LAT ~ loc$LONG, pch=19, col="blue", lwd=2)
dev.off()

#----------------------------------

png(file="NZ2010.png", width=900, height=600)
par(mfrow=c(1,2))

map("nzHires", xlim=c(165,180), ylim=c(-48,-34))
map.axes()
mtext("2010 Earthquakes of magnitude less than 4\n")
points(GNS$LAT[GNS$MAG<4]~GNS$LONG[GNS$MAG<4], pch=19)

map("nzHires", xlim=c(165,180), ylim=c(-48,-34))
map.axes()
mtext("2010 Earthquakes of magnitude greater than 4\n")
points(GNS$LAT[GNS$MAG>4]~GNS$LONG[GNS$MAG>4], pch=19, col="orange")
points(GNS$LAT[GNS$MAG>5]~GNS$LONG[GNS$MAG>5], pch=19, col="blue")
legend(166, -34, legend=c("> 4", "> 5"), pch=19, col=c("orange", "blue"), title="Magnitude")

dev.off()

#----------------------------------

png(file="histograms.png", width=900, height=600)
par(mfrow=c(1,2))
hist(GNS$MAG, main="Magnitude of New Zealand earthqukes \n detected in 2010 to 8 September", xlab="Magnitude (Richter scale)", col="grey25")

hist(GNS$DEPTH, main="Depth of New Zealand earthquakes \ndetected in 2010 to 8 September", xlab="Depth (km)", col="grey25")
dev.off()

plot(cant$MAG ~cant$DEPTH)

plot(cant$DEPTH ~cant$Date, type="n")
points(cant$DEPTH ~cant$Date, pch=21, cex=cant$MAG, bg=cant$MAG)


mod<-glm(MAG[-c(1,2)]~time_diff[-c(1,2)], data=after)


plot(after$MAG ~ after$Date, type="l")
abline(mod$coefficients[1], mod$coefficients[2])