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

#Read in data
GNS<-read.csv("NZ_earthquakes.csv")

#Format Date properly
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")

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


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

#New variables
after$time_diff<-after$Date-after$Date[2]   #Time afterward
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,])

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

plot(GNS$LAT~GNS$LONG, type="n")
points(GNS$LAT~GNS$LONG, cex=GNS$MAG)


plot(GNS$MAG ~ GNS$Date, type="p", pch=18, ylab="Magnitude", xlab="Date")
points(cant$MAG ~ cant$Date, pch=18, col="blue")
legend(1283609257, 7, legend=c("Canterbury", "Elsewherein NZ"), pch=18, col=c("blue", "black"))

plot(cant$MAG~cant$Date, type="l", xlab="Time", ylab="Magnitude")

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

plot(cant$LAT ~ cant$LONG)


#######################
#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="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()