R plotting example: Polar Drift Direction

DriftDirection.R

# Data preparation

# Not included in this code snippet to respect the current proprietary nature of this research.
# Once research is to be reproduced, data preparation procedures will be made available.

#Plotting

par(mar=c(5,4,4,5)+.1)
plot(newyears, pathtangent, type="l", col="white",xlim=c(1890, 2030), ylim=c(30, 225), xaxt="n", yaxt="n", ylab="", xlab="")
axis(1, at=seq(from=1890, to=2030, by=10))
axis(2, at=seq(from=40, to=225, by=20))
abline(v=seq(from=1880, to=2035, by=1), col="green", lwd=0.3)
abline(v=seq(from=1890, to=2030, by=5), col="green")
abline(h=seq(from=20, to=235, by=2), col="green", lwd=0.3)
abline(h=seq(from=30, to=225, by=10), col="green")
par(new=T)
plot(newyears, pathtangent, type="l", col="blue",xlim=c(1890, 2030), ylim=c(30, 225), xaxt="n", yaxt="n", ylab="", xlab="")
par(new=T)
plot(futureyears2[futureyears2<2030.05], futurepathtangent[futureyears2<2030.05], type="l", col="blue", lty="22", xlim=c(1890, 2030), ylim=c(30, 225), xaxt="n", yaxt="n", ylab="", xlab="")

points(smallquakes, placesmallquakes, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=22, cex=0.7, col="red")
points(medquakes, placemedquakes, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=22, cex=1.1, col="red")
points(medquakes, placemedquakes, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=20, cex=1, col="black")
points(bigquakes, placebigquakes, xlim=c(1890, 2010), ylim=c(30, 225), xlab="", ylab="", type="p", pch=22, cex=1.9, col="red")
points(bigquakes, placebigquakes, xlim=c(1890, 2010), ylim=c(30, 225), xlab="", ylab="", type="p", pch=15, cex=1.5, col="black")

points(smallvolcanos, smallplace, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=2, cex=0.7, col="red")
points(medvolcanos, medplace, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=2, cex=1.1, col="red")
points(medvolcanos, medplace, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=20, cex=1, col="black")
points(largevolcanos, largeplace, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=2, cex=1.7, col="red")
points(largevolcanos, largeplace, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=17, cex=1.5, col="black")

text(medquakes[2], placemedquakes[2], labels="(2 ea)", pos=2, cex=0.7)
text(medquakes[11]+1, placemedquakes[11]-5, labels="(2 ea of    )", pos=4, cex=0.7)
par(new=T)
points(medquakes[11]+11.5, placemedquakes[11]-5, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=22, cex=1.1, col="red")
points(medquakes[11]+11.5, placemedquakes[11]-5,xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=20, cex=1, col="black")
text(medquakes[24], placemedquakes[24]-5, labels="(1 ea of    )", pos=4, cex=0.7)
points(medquakes[24]+10.5, placemedquakes[24]-5, xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=22, cex=1.1, col="red")
points(medquakes[24]+10.5, placemedquakes[24]-5,xlim=c(1890, 2030), ylim=c(30, 225), xlab="", ylab="", type="p", pch=20, cex=1, col="black")

lines(c(medquakes[2], medquakes[2]-1.5), c(placemedquakes[2], placemedquakes[2]))
lines(c(medquakes[11], medquakes[11]+2.8), c(placemedquakes[11], placemedquakes[11]-5))
lines(c(medquakes[24], medquakes[24]+1.5), c(placemedquakes[24], placemedquakes[24]-4))

mtext("MGC (E. Longitude)", side=2, line=3)
mtext("YEAR", side=1, line=3)

categories1 <- expression(paste(bold("Traces")), paste("Path of MGCs derived from pole data through Jan 18, 2011:"), paste("Averaging intervals ", t%+-%0.15, " yr pre-1900 ", t%+-%0.05, " yr post-1900"), paste("Path of MGCs forecasted by median of ARIMA model to peak at 2021.3, then mirrored at peak"))
categories2 <- expression(paste(bold("Earthquake Size")), paste("8.0", "" <= "","M",""<= 8.3), paste(8.4<= "","M",""<= 8.7), paste(8.8<= "","M",""<= 9.5), paste(bold("Eruption Size")), paste(3.9<= "","VEI",""<= 4), paste(4< "","VEI",""<= 5), paste(5< "","VEI",""<= 6))

legend(x=1912.5, y=50, categories1, title.col="white", col=c("white", "blue", "white", "blue"), pch=c(-1, -1,-1,-1), lty=c("solid", "solid", "solid", "22"), lwd=c(1, 2,1,2), bty="o", bg="white", cex=0.6)
legend(x=1997, y=50, categories2, title.col="white", col=c("white", "white", "black", "black", "white", "white", "black", "black"), pch=c(-1, -1, 20, 15, -1, -1, 20, 17), pt.cex=c(1,1, 1, 1.5, 1,1,1,1.5), bty="o", bg="white", cex=0.6, ncol=2)
legend(x=1997, y=50, categories2, title.col="white", col=c("white", "red", "red", "red", "white", "red", "red", "red"), pch=c(-1, 22, 22, 22, -1, 2, 2, 2), pt.cex=c(1, 0.7, 1.1, 1.9, 1, 0.7, 1.2, 1.7),bty="o", bg="#FF000000", cex=0.6, ncol=2)


par(new=T)
plot(c(1890, 2030), c(-0.008, 0.008), type="n", xlim=c(1890, 2010), ylim=c(-0.008, 0.008), ylab="", xlab="", yaxt="n", xaxt="n")
rect(1970, 0.0068, 2010, 0.0085, col="white")

text(1970, 0.0081, labels="POLAR DRIFT DIRECTION", pos=4, cex=0.7)
text(1970, 0.0077, labels="In Meridian Great Circles", pos=4, cex=0.7)
text(1970, 0.0073, labels="With Great Earthquakes and Large Volcanic Eruptions", pos=4, cex=0.6)
text(1970, 0.0070, labels="1890 to 2030", pos=4, cex=0.6)
copyright.draw(" IfCG ", z, "BY-ND-2011", .08, .045, .4, fontsize=8)