# Sunrise and sunset times for Chicago IL # Used to illustrate graphics in R # Author: Gene Leynes # Date: February 2012 ############################################################################### #library(geneorama) library(XML) library(RColorBrewer) ## You will have to set up your own directory ## Make sure that "Functions" is a sub directory with the dependent functions setwd('~/Dropbox/R/+Projects/XLConnect') ## Delete all items from the workspace rm(list=ls()) ##----------------------------------------------------------------------------- ## Load data that has been "scraped" from the internet ## Combine tables into a single data frame ##----------------------------------------------------------------------------- load(file='DataFullChicago.RData') str(dat) ls() ##----------------------------------------------------------------------------- ## Convert time data to a "fraction of the day" type of data ##----------------------------------------------------------------------------- ## I want to treat "time of day" as a value, so that I can view it over days. ## As this plot reveals the "time" data isn't handled that way by default: plot(dat[,1:2]) ## So this function will convert it to a value ## For example, noon will become .5, and 3 PM will be 15/24 (or 0.625) source('Functions/ConvertDateTimeToDecimal.R') TimeColumns = grep('POSIX', sapply(dat, class)) dat[,TimeColumns] = sapply(dat[,TimeColumns], ConvertDateTimeToDecimal) str(dat) ## Now this plot makes a little more sense: plot(dat[,1:2]) ##----------------------------------------------------------------------------- ## Fix the daylight savings time ##----------------------------------------------------------------------------- ## The last plot shows a big disconnect when daylight savings time happens ## To fix that I'm going to use a function that I wrote source('Functions/FixDaylightSavings.R') ## Test the function: par(mfrow=c(1,2)) plot(dat$Solar_Noon_Time, main='Solar Noon (not adjusted)') plot(FixDaylightSavings(dat$Solar_Noon_Time), main='Solar Noon (adjusted)') par(mfrow=c(1,1)) ## Apply function to a copy of "dat", on all time date columns datAdj = dat datAdj[,TimeColumns] = sapply(datAdj[,TimeColumns], FixDaylightSavings) ##----------------------------------------------------------------------------- ## Column Names (for reference) ##----------------------------------------------------------------------------- ## clipper is in the "geneorama" package, which is not yet published ## Here is the definition of clipper # clipper=function(x){ # write.table(x,"clipboard",sep="\t",col.names=TRUE,row.names=FALSE) #} # clipper(colnames(dat)) #Astronomical_Twilight_Starts #Astronomical_Twilight_Ends #Nautical_Twilight_Starts #Nautical_Twilight_Ends #Civil_Twilight_Starts #Civil_Twilight_Ends #Sunrise #Sunset #Azimuth_Sunrise #Azimuth_Sunset #Day_Length #Day_Length_Difference #Solar_Noon_Time #Solar_Noon_Altitude #Distance_10e6_KM ##----------------------------------------------------------------------------- ## Initial Plots ##----------------------------------------------------------------------------- plot(dat$Astronomical_Twilight_Starts, main='Twilight Starts') plot(dat$Astronomical_Twilight_Ends, main='Twilight Ends') matplot(dat[,c(2,3)], type='l', main='Twilight Start and End (Astronomical)') matplot(dat[,c(4,5)], type='l', main='Twilight Start and End (Nautical)') matplot(dat[,c(4,5)], type='l', main='Twilight Start and End (Nautical)') matplot(dat[,c(2,3,4,5)], type='l', main='Twilight Start and End (Astro and Naut)') matplot(dat[,c(2,3,4,5)], type='l', main=paste('Twilight Start and End', '(Astronomical and Nautical)',sep='\n')) matplot(dat[,c(2,3,4,5,6,7)], type='l', main=paste('Twilight Start and End', '(Astronomical, Nautical, and Civil)',sep='\n')) ##----------------------------------------------------------------------------- ## Looking at color palates ##----------------------------------------------------------------------------- #plot100colors= function (x) { # dat = data.frame(x = rep(1:10, each = 10), y = rep(1:10, 10)) # dat$y = 11 - dat$y # z = matrix(x:(x + 99), 10, 10, byrow = TRUE) # plot(dat, cex = 5, col = colors()[z], pch = 16, yaxt = "n", # xaxt = "n", xlab = "", ylab = "") # points(dat, cex = 5, pch = 21, lwd = 1) # text(dat$x, dat$y - 0.3, colors()[z], cex = 0.6) # text(dat$x, dat$y - 0.4, z, cex = 0.6) #} #plot100colors(1) #plot100colors(101) #plot100colors(201) #plot100colors(301) #plot100colors(401) #plot100colors(501) #plot100colors(601) #barplot(rep(1,5), col=brewer.pal(9,'YlGnBu')) ##----------------------------------------------------------------------------- ## Plotting the extent of the day (with daylight savings) ##----------------------------------------------------------------------------- cols=brewer.pal(9,'YlGnBu') #barplot(rep(1,length(cols)), col=cols) plot(dat[,2], type='n', ylim=c(0,1), ann=F, axes=F) xx = 1:nrow(dat) polygon(rbind(c(0,0), c(0,1), c(nrow(dat),1), c(nrow(dat),0)), col=cols[9], border = NA) polygon(rbind(cbind(xx,dat[,2]),cbind(rev(xx),rev(dat[,3]))), col=cols[7], border = NA) polygon(rbind(cbind(xx,dat[,4]),cbind(rev(xx),rev(dat[,5]))), col=cols[5], border = NA) polygon(rbind(cbind(xx,dat[,6]),cbind(rev(xx),rev(dat[,7]))), col=cols[4], border = NA) polygon(rbind(cbind(xx,dat[,8]),cbind(rev(xx),rev(dat[,9]))), col=cols[1], border = NA) sapply(seq(3,21,3)/24, function(x) lines(xx,rep(x,nrow(dat)),col='grey75', lty=1)) lines(xx,dat$Solar_Noon_Time,col='red', lty=2, lwd=2) text(300, .57, 'Actual solar noon', col='red') yvals = seq(0,24,3)/24 ylabs = c('12:00 AM','3:00 AM','6:00 AM','9:00 AM', '12:00 PM','3:00 PM','6:00 PM','9:00 PM','12:00 PM') axis(2, at=yvals, labels=ylabs, col.axis="grey25", tick=F, las=2, cex.axis=.7) months = c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec') axis(1, at=15+seq(0,330,30), labels=months, tick=T, las=1, cex.axis=.7) #plot(dat[,2], type='n', ylim=c(0,1), ann=F, axes=F) mtext('Daylight with daylight savings', side=3, line=2.5, cex=1.5) mtext('Gradients indicate twilight (astronomical, nautical, and civil)', side=3, line=1, cex=1) mtext('source: http://www.timeanddate.com/', side=1, line=4, cex=1) mtext(paste('Times shown are for ', attr(dat,'location')), side=1, line=3, cex=1) ##----------------------------------------------------------------------------- ## Plotting the extent of the day (without daylight savings) ##----------------------------------------------------------------------------- cols=brewer.pal(9,'YlGnBu') plot(datAdj[,2], type='n', ylim=c(0,1), ann=F, axes=F) xx = 1:nrow(datAdj) polygon(rbind(c(0,0), c(0,1), c(nrow(datAdj),1), c(nrow(datAdj),0)), col=cols[9], border = NA) polygon(rbind(cbind(xx,datAdj[,2]),cbind(rev(xx),rev(datAdj[,3]))), col=cols[7], border = NA) polygon(rbind(cbind(xx,datAdj[,4]),cbind(rev(xx),rev(datAdj[,5]))), col=cols[5], border = NA) polygon(rbind(cbind(xx,datAdj[,6]),cbind(rev(xx),rev(datAdj[,7]))), col=cols[4], border = NA) polygon(rbind(cbind(xx,datAdj[,8]),cbind(rev(xx),rev(datAdj[,9]))), col=cols[1], border = NA) sapply(seq(3,21,3)/24, function(x) lines(xx,rep(x,nrow(datAdj)),col='grey75', lty=1)) lines(xx,datAdj$Solar_Noon_Time,col='red', lty=2, lwd=2) text(300, .57, 'Actual solar noon', col='red') yvals = seq(0,24,3)/24 ylabs = c('12:00 AM','3:00 AM','6:00 AM','9:00 AM', '12:00 PM','3:00 PM','6:00 PM','9:00 PM','12:00 PM') axis(2, at=yvals, labels=ylabs, col.axis="grey25", tick=F, las=2, cex.axis=.7) months = c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec') axis(1, at=15+seq(0,330,30), labels=months, tick=T, las=1, cex.axis=.7) #plot(datAdj[,2], type='n', ylim=c(0,1), ann=F, axes=F) mtext('Daylight *without* daylight savings', side=3, line=2.5, cex=1.5) mtext('Gradients indicate twilight (astronomical, nautical, and civil)', side=3, line=1, cex=1) mtext('source: http://www.timeanddate.com/', side=1, line=4, cex=1) mtext(paste('Times shown are for ', attr(dat,'location')), side=1, line=3, cex=1) #Nautical_Twilight_Starts #Nautical_Twilight_Ends #Civil_Twilight_Starts #Civil_Twilight_Ends #Sunrise #Sunset #Azimuth_Sunrise #Azimuth_Sunset #Day_Length #Day_Length_Difference #Solar_Noon_Time #Solar_Noon_Altitude #Distance_10e6_KM plot(dat$Distance_10e6_KM, main='Distance to sun') plot(dat$Solar_Noon_Altitude, main='Angle of sun at solar noon') plot(dat$Day_Length, main='Day Length') plot(dat$Day_Length_Difference, main='Difference in day lengths') ##----------------------------------------------------------------------------- ## Simple plot of change in length of day ##----------------------------------------------------------------------------- plot(dat$Day_Length_Difference * 24 * 60, axes=F, ann=F, pch=16, cex=.7) abline(h=0) months = c('Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec') axis(1, at=15+seq(0,330,30), labels=months, tick=T, las=1, cex.axis=.7) axis(2) mtext('Change in length of day (in minutes)', side=3, line=2, cex=1.5) mtext('source: http://www.timeanddate.com/', side=1, line=4, cex=1) mtext(paste('Times shown are for ', attr(dat,'location')), side=1, line=3, cex=1)