[This article was first published on DataDebrief, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 their blog: DataDebrief.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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