Plotting time vs date in R

October 2, 2010
By

(This article was first published on DataDebrief, and kindly contributed to R-bloggers)

Having done the plot with python+matplotlib, thought I would have a go reproducing it in R, using only the builtin "plot". So, just to recap - this is a plot of the sun times (rise/set,twilight and blinding) as generated by the great python library pyephem. The R code reads in a csv file as output from a modified version of the python code used in my original post.

For completeness, the csv generation code is below:

import ephem
import datetime, math
import pylab

place = ephem.city('Melbourne')

start_date = datetime.datetime(2009,12,1,12)
end_date = datetime.datetime(2011, 1, 31,12)

base_offset = '0'
twilight_offset = '-6:00:0.0' # "twilight" = centre of the sun is -6deg ideal horizon
eyeline_offset = '15:34:0.0' # arbitrary +15deg

sun = ephem.Sun(place)

dates = []
sunrise = []
sunset = []
firstlight = []
lastlight = []
firsteyel = []
lasteyel = []

numdays = (end_date - start_date).days
dates = [start_date + datetime.timedelta(days=i) for i in xrange(numdays+1)]
dates.sort()

def dt2m(dt):
return (dt.hour*60) + dt.minute

def m2hm(x):
h = int(x/60)
m = int(x%60)
return '%(h)02d%(m)02d' % {'h':h,'m':m}

sunrise = map(lambda x:dt2m(ephem.localtime(place.next_rising(sun,start=x))),dates)
sunset = map(lambda x:dt2m(ephem.localtime(place.next_setting(sun,start=x))),dates)

place.horizon = twilight_offset
firstlight = map(lambda x:dt2m(ephem.localtime(place.next_rising(sun,start=x))),dates)
lastlight = map(lambda x:dt2m(ephem.localtime(place.next_setting(sun,start=x))),dates)

place.horizon = eyeline_offset
firsteyel = map(lambda x:dt2m(ephem.localtime(place.next_rising(sun,start=x))),dates)
lasteyel = map(lambda x:dt2m(ephem.localtime(place.next_setting(sun,start=x))),dates)

writer = open("suntimes.csv", "w")
writer.write("date,firstlight,sunrise,firsteyel,lasteyel,sunset,lastlight\n")
for n in xrange(numdays):
writer.write(str(dates[n]) +","+ (m2hm(firstlight[n]))+"," \
+ m2hm(sunrise[n])+"," + m2hm(firsteyel[n])+"," + m2hm(lasteyel[n])+"," \
+ m2hm(sunset[n])+"," + m2hm(lastlight[n])+"\n")
writer.close()


I chose to leave the format of the times, as they resemble the format of the Geoscience Australia times.

This is how the csv file looks:
date,firstlight,sunrise,firsteyel,lasteyel,sunset,lastlight
2009-12-01 00:00:00,0518,0551,0719,1858,2026,2059
2009-12-02 00:00:00,0518,0551,0719,1859,2027,2100
2009-12-03 00:00:00,0518,0551,0719,1859,2028,2101
2009-12-04 00:00:00,0518,0551,0719,1900,2029,2102
2009-12-05 00:00:00,0518,0550,0719,1901,2030,2103


And here's the resulting graph:


Here's the R code:

suntimes <- read.csv("suntimes.csv")

suntimes$date <- sub(' 00:00:00$', '', suntimes$date) #strip off trailing 0
suntimes$date <- as.Date(suntimes$date,"%Y-%m-%d") # make real date so we can axis.date
suntimes$firstlight <- sub('^([[:digit:]]{3})$', '0\\1', suntimes$firstlight) #pad 0
suntimes$sunrise <- sub('^([[:digit:]]{3})$', '0\\1', suntimes$sunrise)
suntimes$firsteyel <- sub('^([[:digit:]]{3})$', '0\\1', suntimes$firsteyel)
suntimes$lasteyel <- sub('^([[:digit:]]{3})$', '0\\1', suntimes$lasteyel)
suntimes$sunset <- sub('^([[:digit:]]{3})$', '0\\1', suntimes$sunset)
suntimes$lastlight <- sub('^([[:digit:]]{3})$', '0\\1', suntimes$lastlight)

#calc as minutes from midnight
suntimes$firstlighttimes <- as.integer(substr(suntimes$firstlight,1,2))*60 + as.integer(substr(suntimes$firstlight,3,4))
suntimes$sunrisetimes <- as.integer(substr(suntimes$sunrise,1,2))*60 + as.integer(substr(suntimes$sunrise,3,4))
suntimes$firsteyeltimes <- as.integer(substr(suntimes$firsteyel,1,2))*60 + as.integer(substr(suntimes$firsteyel,3,4))
suntimes$lasteyeltimes <- as.integer(substr(suntimes$lasteyel,1,2))*60 + as.integer(substr(suntimes$lasteyel,3,4))
suntimes$sunsettimes <- as.integer(substr(suntimes$sunset,1,2))*60 + as.integer(substr(suntimes$sunset,3,4))
suntimes$lastlighttimes <- as.integer(substr(suntimes$lastlight,1,2))*60 + as.integer(substr(suntimes$lastlight,3,4))

plot(suntimes$date,suntimes$firstlighttimes,type="l", lwd=1,col='blue',ylim=c(240,1330), axes= F, xlab= "Dates", ylab= "Time")
# make grid first, then overlay points
abline(h=(seq(from=240,to=1330, by=20)), lwd =0.1, lty="dotted", col='#fafafa')
abline(h=(seq(from=240,to=1330, by=60)), lwd =0.1, lty="dotted", col='#fceae4')
abline(v=(seq(from=min(suntimes$date),to=max(suntimes$date),"month")), lwd =0.1, lty="dotted", col='#fceae4')

points(suntimes$date,suntimes$firstlighttimes,type='l', lwd=1,col='blue')
#points(suntimes$date,suntimes$sunrisetimes,type='p',pch=20, cex = 0.75, lwd=0.5,col='#FF0000')
points(suntimes$date,suntimes$sunrisetimes,type='l', lwd=2,col='#FF0000')
points(suntimes$date,suntimes$firsteyeltimes,type='l', lwd=1,col='orange')
points(suntimes$date,suntimes$lasteyeltimes,type='l', lwd=1,col='orange')
#points(suntimes$date,suntimes$sunsettimes,type='p',pch=20, cex = 0.75, lwd=0.75, col='#FF0000')
points(suntimes$date,suntimes$sunsettimes,type='l', lwd=2, col='#FF0000')
points(suntimes$date,suntimes$lastlighttimes,type='l', lwd=1,col='blue')

#customise X Axis
axis.Date(side=1, at=seq(from=min(suntimes$date),to=max(suntimes$date),"days"),col.ticks='gray', format="%b", labels="")
axis.Date(side=1, at=seq(from=min(suntimes$date),to=max(suntimes$date),"month"), col.ticks='red', format="%b\n%Y" )

#customise Y Axis
yaxisMajorTicks.hours <- seq(from=240,to=1330, by=60) # whole hours from 4am to 10pm
yaxisMajorTicks.hournames <- as.character(yaxisMajorTicks.hours/60) #eg 240/60=4am
yaxisMajorTicks.hournames <- sub('^(.+)$', '\\1:00', yaxisMajorTicks.hournames) #append ":00"

yaxisMinorTicks.tens <- seq(from=240,to=1330, by=10)
yaxisMinorTicks.tensnames <- as.character(yaxisMinorTicks.tens)
yaxisMinorTicks.tensnames <- sub('^(.+)$', "", yaxisMinorTicks.tensnames)
# draw major over minor ticks
axis(side=2, at=yaxisMinorTicks.tens, col.ticks='gray',labels=yaxisMinorTicks.tensnames)
axis(side=2, at=yaxisMajorTicks.hours, col.ticks='red',labels=yaxisMajorTicks.hournames, las=2)
box()
#legend max(yaxisMajorTicks.hours),c(1,2,1,2,1)pch=c(NA,20,NA,20,NA)
legend("right","center", c("firstlight","sunrise","in eyes","sunset","lastlight"), cex=0.75,
col=c("blue","red","orange","red","blue"), lty=1, lwd=2, bty="o", bg='#ffffff');


To leave a comment for the author, please follow the link and comment on his blog: DataDebrief.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.