# Sunrise and sunset times for Chicago IL # Used to illustrate XLConnect package, and graphics in R # Author: Gene ############################################################################### #library(geneorama) #library(XLConnect) 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/XLConnect') rm(list=ls()) ##----------------------------------------------------------------------------- ## Define initial URL and use it to create a set of URLs (one for each month) ##----------------------------------------------------------------------------- ## Chicago sunriseURLIL = paste( 'http://www.timeanddate.com/worldclock/astronomy.html', '?n=64&month=1&year=2012&obj=sun&afl=-11&day=1', sep='') ## Anchorage sunriseURLAK = paste( 'http://www.timeanddate.com/worldclock/astronomy.html', '?n=18&month=1&year=2012&obj=sun&afl=-11&day=1', sep='') ## Miami sunriseURLFL = paste( 'http://www.timeanddate.com/worldclock/astronomy.html', '?n=156&month=1&year=2012&obj=sun&afl=-11&day=1', sep='') sunriseURLsIL = sapply(1:12, function(x){ NewText = paste('month=',x,sep='') gsub('month=1',NewText,sunriseURLIL) }) sunriseURLsAK = sapply(1:12, function(x){ NewText = paste('month=',x,sep='') gsub('month=1',NewText,sunriseURLAK) }) sunriseURLsFL = sapply(1:12, function(x){ NewText = paste('month=',x,sep='') gsub('month=1',NewText,sunriseURLFL) }) sunriseURLAK ## Original URL sunriseURLsAK ## New URLs based on original sunriseURLsIL ## New URLs based on original sunriseURLsFL ## New URLs based on original ##----------------------------------------------------------------------------- ## Use XML library to get HTML tables ##----------------------------------------------------------------------------- RawHtmlDataIL = lapply(sunriseURLsIL, function(x) readHTMLTable(x, header=T, which=1, stringsAsFactors=F)) RawHtmlDataAK = lapply(sunriseURLsAK, function(x) readHTMLTable(x, header=T, which=1, stringsAsFactors=F)) RawHtmlDataFL = lapply(sunriseURLsFL, function(x) readHTMLTable(x, header=T, which=1, stringsAsFactors=F)) head(RawHtmlDataAK[[1]]) head(RawHtmlDataIL[[1]]) head(RawHtmlDataFL[[1]]) str(RawHtmlDataAK[[1]]) str(RawHtmlDataIL[[1]]) str(RawHtmlDataFL[[1]]) ## A little backup so that I don't have to re-run the code ## I normally run this part manually and commenting out the lines above if(FALSE){ #save(RawHtmlDataIL, file='RawHtmlDataIL.RData') #save(RawHtmlDataAK, file='RawHtmlDataAK.RData') #save(RawHtmlDataFL, file='RawHtmlDataFL.RData') RawHtmlDataIL = loader('RawHtmlDataIL') RawHtmlDataAK = loader('RawHtmlDataAK') RawHtmlDataFL = loader('RawHtmlDataFL') } ##----------------------------------------------------------------------------- ## Use custom function to clean up tables ##----------------------------------------------------------------------------- source('Functions/CleanUpTable_simple.R') CleanHtmlDataIL = lapply(RawHtmlDataIL, CleanUpTable_simple) CleanHtmlDataAK = lapply(RawHtmlDataAK, CleanUpTable_simple) CleanHtmlDataFL = lapply(RawHtmlDataFL, CleanUpTable_simple) str(CleanHtmlDataIL) str(CleanHtmlDataAK) str(CleanHtmlDataFL) ##----------------------------------------------------------------------------- ## Combine tables into a single data frame ##----------------------------------------------------------------------------- datIL = do.call(rbind, CleanHtmlDataIL) datAK = do.call(rbind, CleanHtmlDataAK) datFL = do.call(rbind, CleanHtmlDataFL) str(datIL) str(datAK) str(datFL) ##----------------------------------------------------------------------------- ## Convert time data to a "fraction of the day" type of data ##----------------------------------------------------------------------------- source('Functions/ConvertDateTimeToDecimal.R') TimeColumns = grep('POSIX', sapply(datIL, class)) datIL[,TimeColumns] = sapply(datIL[,TimeColumns], ConvertDateTimeToDecimal) datAK[,TimeColumns] = sapply(datAK[,TimeColumns], ConvertDateTimeToDecimal) datFL[,TimeColumns] = sapply(datFL[,TimeColumns], ConvertDateTimeToDecimal) str(datIL) str(datAK) str(datFL) ##----------------------------------------------------------------------------- ## Initial Plots ##----------------------------------------------------------------------------- colnames(datIL) matplot(1:nrow(datIL), datIL[,c('Sunrise', 'Sunset')], main='Sunrise and Sunset', type='l', col='black', lty=1, ylim=c(0,1)) lines(1:nrow(datAK), datAK[,'Sunrise'], col='blue') lines(1:nrow(datAK), datAK[,'Sunset'], col='blue') lines(1:nrow(datFL), datFL[,'Sunrise'], col='red') lines(1:nrow(datFL), datFL[,'Sunset'], col='red') ##----------------------------------------------------------------------------- ## Plotting the extent of the day (Chicago vs Anchorage) ##----------------------------------------------------------------------------- cols=brewer.pal(9,'YlGnBu') plot(datIL[,2], type='n', ylim=c(0,1), ann=F, axes=F) xx = 1:nrow(datIL) polygon(rbind(c(0,0), c(0,1), c(nrow(datIL),1), c(nrow(datIL),0)), col=cols[5], border = NA) polygon(rbind(cbind(xx,datIL$Sunrise),cbind(rev(xx),rev(datIL$Sunset))), col=cols[1], border = NA) sapply(seq(3,21,3)/24, function(x) lines(xx,rep(x,nrow(datIL)),col='grey75', lty=1)) lines(xx,datIL$Solar_Noon_Time,col='red', lty=2, lwd=2) lines(xx, datAK[,'Sunrise'], col='blue', lwd=2) lines(xx, datAK[,'Sunset'], col='blue', lwd=2) lines(xx,datAK$Solar_Noon_Time,col='blue', lty=2, lwd=2) 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('Chicago Daylight vs. Anchorage Daylight', side=3, line=2.5, cex=1.5) mtext('Anchorage daylight shown in blue', side=3, line=1, cex=1) text(300, .45, 'Chicago solar noon', col='red') text(300, .61, 'Anchorage solar noon', col='blue') mtext('source: http://www.timeanddate.com/', side=1, line=4, cex=1) ##----------------------------------------------------------------------------- ## Plotting the extent of the day (Chicago vs Miami) ##----------------------------------------------------------------------------- cols=brewer.pal(9,'YlGnBu') plot(datIL[,2], type='n', ylim=c(0,1), ann=F, axes=F) xx = 1:nrow(datIL) polygon(rbind(c(0,0), c(0,1), c(nrow(datIL),1), c(nrow(datIL),0)), col=cols[5], border = NA) polygon(rbind(cbind(xx,datIL$Sunrise),cbind(rev(xx),rev(datIL$Sunset))), col=cols[1], border = NA) sapply(seq(3,21,3)/24, function(x) lines(xx,rep(x,nrow(datIL)),col='grey75', lty=1)) lines(xx,datIL$Solar_Noon_Time,col='red', lty=2, lwd=2) lines(xx, datFL[,'Sunrise'], col='orange', lwd=2) lines(xx, datFL[,'Sunset'], col='orange', lwd=2) lines(xx,datFL$Solar_Noon_Time,col='orange', lty=2, lwd=2) 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('Chicago Daylight vs. Miami Daylight', side=3, line=2.5, cex=1.5) mtext('Miami daylight shown in orange', side=3, line=1, cex=1) text(300, .45, 'Chicago solar noon', col='red') text(300, .58, 'Miami solar noon', col='orange') mtext('source: http://www.timeanddate.com/', side=1, line=4, cex=1) plot(datIL$Distance_10e6_KM) plot(datIL$Solar_Noon_Altitude) plot(datIL$Day_Length) plot(datIL$Day_Length_Difference) ##----------------------------------------------------------------------------- ## Simple plot of change in length of day ##----------------------------------------------------------------------------- plot(datIL$Day_Length_Difference * 24 * 60, axes=F, ann=F, pch=16, cex=.7, ylim=c(-7.5, 7.5)) points(datAK$Day_Length_Difference * 24 * 60, col='blue', pch=16, cex=.7) points(datFL$Day_Length_Difference * 24 * 60, col='orange', 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) legend(250, 7.5, c("Anchorage", "Chicago", "Miami"), col = c('blue', 'black', 'orange'), text.col = "black", lty = c(1,1,1), pch = NA, lwd = 2, merge = TRUE, bg = 'gray90') mtext('source: http://www.timeanddate.com/', side=1, line=4, cex=1)