# Data from Land Information New Zealand # See Moon/provenance.txt and Moon/moon.R lowTides <- c("00:55", "01:47", "02:38", "03:29", "04:20", "05:11", "06:05", "07:01", "08:00", "09:01", "10:01", "10:58", "11:51", "00:16", "01:03", "01:44", "02:23", "03:00", "03:36", "04:12", "04:49", "05:29", "06:12", "07:02", "07:59", "09:02", "10:06", "11:08", "12:08", "00:34", "01:29") phases <- data.frame(date=c("2010-01-01 08:13", "2010-01-07 23:40", "2010-01-15 08:12", "2010-01-23 23:54", "2010-01-30 07:18"), phase=c("Full", "3Q", "New", "1Q", "Full")) phases$date <- as.POSIXct(as.character(phases$date)) lowTideDate <- as.POSIXct(paste("2010-01-", sprintf("%02d", 1:31), " ", lowTides, sep="")) lowTideHour <- as.POSIXct(paste("2010-01-01 ", lowTides, sep="")) mainHours <- ISOdatetime(2010, 1, 1, c(0, 4, 8, 12), rep(0, 4), rep(0, 4)) library(jpeg) moon <- jpeg::readJPEG(system.file("extra", "GPN-2000-000473.jpg", package="RGraphics")) moonPhase <- function(x, y, phase, size=.05) { # size is in inches n <- 17 angle <- seq(0, 2*pi, length=n) xx <- x + cos(angle)*xinch(size) yy <- y + sin(angle)*yinch(size) if (phase == "New") fill <- "black" else fill <- "white" polygon(xx, yy, col=fill) if (phase == "1Q") polygon(xx[(n/4):(n*3/4) + 1], yy[(n/4):(n*3/4) + 1], col="black") if (phase == "3Q") polygon(xx[c(1:(n/4 + 1), (n*3/4 + 1):n)], yy[c(1:(n/4 + 1), (n*3/4 + 1):n)], col="black") } # Original image from NASA # http://grin.hq.nasa.gov/ABSTRACTS/GPN-2000-000473.html rasterMoon <- jpeg::readJPEG("Moon/GPN-2000-000473.jpg") par(pin=c(3.5, 1.75), oma=c(0, 3, 0, 0), xaxs="i", yaxs="i", cex=.7) plot.new() rect(0, 0, 1, 1, col="black") rasterImage(rasterMoon, .25, 0, .75, 1*813/703) par(new=TRUE, xaxs="r", yaxs="r", las=1) plot(lowTideDate, lowTideHour, type="n", ylim=range(mainHours), axes=FALSE, ann=FALSE) # dashed reference lines abline(v=phases$date, col="white", lty="dashed") for (subset in list(1:13, 14:29, 30:31)) { lines(lowTideDate[subset], lowTideHour[subset], lwd=2, col="white") points(lowTideDate[subset], lowTideHour[subset], pch=16, col="white") } box() axis.POSIXct(1, lowTideDate) axis.POSIXct(2, at=mainHours, format="%H:%M") mtext("Time of Low Tide (NZDT)", side=2, line=4, las=0, cex=.7) mtext("Auckland, New Zealand January 2010", side=1, line=3, cex=.7) axis(3, at=phases$date, labels=FALSE) par(xpd=NA) ymax <- par("usr")[4] for (i in 1:nrow(phases)) moonPhase(phases$date[i], ymax + yinch(.2), phases$phase[i]) mtext("Phases of the Moon", side=3, line=3, cex=.7)