Plotting time vs date in R

October 2, 2010

(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 ='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)]

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")
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")

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:

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
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)
#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 their blog: DataDebrief. offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, 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.


Mango solutions

plotly webpage

dominolab webpage

Zero Inflated Models and Generalized Linear Mixed Models with R

Quantide: statistical consulting and training




CRC R books series

Six Sigma Online Training

Contact us if you wish to help support R-bloggers, and place your banner here.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)