Plotting time vs date in python and matplotlib

Having to commute to the city each day and frequently getting the blinding sun, I thought as a fun exercise, I could plot the sunrise/sets and thus estimate with some accuracy when the sun is most blinding (using a purely subjective (which seems to fit my anecdotal observations) value for the degrees above the horizon).

The essentials (using python 2.6) are pyephem (astronomical library), numpy (scientific library), and matplotlib (2D graphics and charting library).

Calculating the times is very easy using pyephem, set your observer object to the place you want and off you go, calling next_rising or next_setting, they are returned as datetime objects when wrapped in ephem's "localtime" function (otherwise they're custom dates).

The hard part? plotting times on the Y axis using matplotlib ... if you do it with datetime objects. I gave up after trying a few different ways, and resorted to a very easy workaround.

Plot the times as integer minutes : 0 being midnight, 60 = 1am ... which is 24hr * 60. Then format the Y axis as Hour:Minute using a custom formatter:

from matplotlib.ticker import FuncFormatter as ff


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

ax.yaxis.set_major_formatter(ff(m2hm))

To get sensible tick marks, I set the major ticks to an hour

ax.yaxis.set_major_locator(pylab.MultipleLocator(60))

The minor are set to 15 minutes though I left off showing the labels for the Y minor as it clutters the graph, but left it commented in the code.

import ephem
import datetime, math
import pylab
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from matplotlib.ticker import FuncFormatter as ff

place = ephem.city('Melbourne')

start_date = datetime.datetime(2009,12,1,12) # yyyy,m,d,h - midday is best
end_date = datetime.datetime(2011, 1, 31,12) # midday Jan 31st, 2011

twilight_offset = '-6:00:0.0' # "twilight" = centre of the sun is -6deg ideal horizon
eyeline_offset = '15:00: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, i):
    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)

fig = plt.figure()

ax = fig.add_subplot(111)
ax.plot_date(dates, firstlight, ':',label=' First light (' + twilight_offset + ')', color='#06544e')
ax.plot_date(dates, sunrise, linestyle='-',label=' Sunrise', color='red')
ax.plot_date(dates, firsteyel, linestyle='-', label=' In eyes ('+ eyeline_offset + ')', color='orange')
ax.plot_date(dates, lasteyel, linestyle='-', color='orange')
ax.plot_date(dates, sunset, linestyle='-', label=' Sunset', color='red')
ax.plot_date(dates, lastlight, ':',label=' Last light (' + twilight_offset + ')', color='blue')

ax.set_xlim([dates[0],dates[len(dates)-1]])
ax.xaxis.set_major_locator(mdates.AutoDateLocator())
ax.xaxis.set_minor_locator(mdates.DayLocator())
#ax.xaxis.grid(False)
ax.xaxis.grid(color='#ffcf06', linestyle=':', linewidth=0.5)

ax.set_ylim(240,1330) #4am to 10pm
ax.yaxis.set_major_locator(pylab.MultipleLocator(60))
ax.yaxis.set_minor_locator(pylab.MultipleLocator(15))
ax.yaxis.set_major_formatter(ff(m2hm))
#ax.yaxis.set_minor_formatter(ff(m2hm))
ax.yaxis.grid(color='r', linestyle=':', linewidth=0.5)

labels = ax.get_xticklabels()
for label in labels:
    label.set_rotation(30)
plt.legend(loc='best')
plt.xlabel('Date' )
plt.ylabel('Hour (24hr)' )
plt.title("Sun times\n" + str(start_date) + " - " + str(end_date) + "\n (localised for " + place.name + ')')
plt.show()