Moon phases with PyEphem

Calculating all sorts of astronomical data including the phases of the moon is a piece of cake if you use the right tools. Combining the fantastic PyEphem package with the powers of Gnuplot and a suitable font will provide you with an actractive illustration of the path of the moon through the sky.

Let us first have a look at what we may achieve. The illustration shows the path of the moon through the sky as calculated for the Netherlands on January 11, 2011. The path of the sun is displayed as well for comparison. The horizontal axis shows the compass direction while the vertical axis shows the elevation above the horizon. The path of the moon is plotted with a resolution of 15 minutes and only for those parts that are above the horizon. On top of the path we have plotted an image indicating its phase for every hour together with the time. Now how did we put this together?

The essential ingredients are the PyEphem package,Gnuplot and a suitable font. The first thing we have to do is calculate the position of the moon as seen from our location and its phase. This is actually quite straightforward:

from datetime import date
from math import radians as rad,degrees as deg

import ephem

g = ephem.Observer()
g.name='Somewhere'
g.lat=rad(52.0)  # lat/long in decimal degrees
g.long=rad(5.2)

m = ephem.Moon()

g.date = date.today()# local time zone, I'm in UTC+1
g.date -= ephem.hour # always everything in UTC

for i in range(24*4): # compute position for every 15 minutes
    m.compute(g)

    nnm = ephem.next_new_moon(g.date)
    pnm = ephem.previous_new_moon(g.date)
    # for use w. moon_phases.ttf A -> just past  newmoon,
    # Z just before newmoon
    # '0' is full, '1' is new
    # note that we cannot use m.phase as this is the percentage of the moon
    # that is illuminated which is not the same as the phase!
    lunation=(g.date-pnm)/(nnm-pnm)
    symbol=lunation*26
    if symbol < 0.2 or symbol > 25.8 :
        symbol = '1'  # new moon
    else:
        symbol = chr(ord('A')+int(symbol+0.5)-1)

    print(ephem.localtime(g.date).time(), deg(m.alt),deg(m.az),
      ephem.localtime(g.date).time().strftime("%H%M"),
      m.phase,symbol)
    g.date += ephem.minute*15

If you run the script a list of numbers is produced that may be processed by Gnuplot. The list looks something like:

00:00:00.000002 7.91511093193 276.570014287 0000 45.1640434265 G
00:15:00.000002 5.7418040753 279.412684378 0015 45.262512207 G
00:29:59.000002 3.61571815838 282.256638545 0029 45.3610076904 G
00:44:59.000002 1.57919683222 285.110182297 0044 45.4595336914 G
...

The first column is the time, the second the elevation above the horizon, the third the azimuth (the compass direction), the forth the time again in a compact HHMM format and the final column holds a letter that will be used to display the phase of the moon with the help of the moon_phases.ttf font.

Assuming you saved the data in a file m.dat the following bit of gnuplot will create a picture called moon.png similar to the one at the top of the page (if you saved it as moon.plot you can run it with gnuplot moon.plot)

set terminal png truecolor nocrop font "/usr/share/fonts/truetype/ttf-dejavu/DejaVuSans.ttf" 12 size 640,480
set output 'moon.png'
set grid xtics nomxtics ytics nomytics noztics nomztics \
 nox2tics nomx2tics noy2tics nomy2tics nocbtics nomcbtics
set grid layerdefault   linetype 0 linewidth 1.000,  linetype 0 linewidth 1.000
set xtics border in scale 1,0.5 mirror norotate  offset character 0, 0, 0
set xtics 45 norangelimit
set ytics border in scale 1,0.5 mirror norotate  offset character 0, 0, 0
set ytics autofreq  norangelimit
set xlabel "Azimuth (compass direction, in degrees)"
set xlabel  offset character 0, 0, 0 font "" textcolor lt -1 norotate
set xrange [ 45.0000 : 315.000 ] noreverse nowriteback
set ylabel "Elevation (degrees above horizon)"
set ylabel  offset character 0, 0, 0 font "" textcolor lt -1 rotate by -270
set yrange [ 0.000000 : 90.0000 ] noreverse nowriteback
plot '/tmp/m.dat'u 3:2 w l lc 4 lw 2 t "moon path", '/tmp/m.dat' u 3:2:(10) every 4 w circles t "" lc rgb "#ddddff", '/tmp/m.dat' u 3:2:6 every 4 w labels offset character 0,character -1 font "/home/michel/sitescripts/sunrise/moon_phases.ttf,30"   t "", '/tmp/m.dat'u 3:2:4 every 4 w labels font "/usr/share/fonts/truetype/ttf-dejavu/DejaVuSans.ttf,8" tc rgb "#dd0000" t ""
#    EOF

The crucial bit is the final (and very long) plot. It draws several subplots:

• a line representing the path of the moon
• a series of circles with a radius of 10 in pale blue to be used as a backdrop for the moon phase characters
• a series of labels that represent the phases of the moon (you must explicitely point to the location of the font)
• and a final set of labels that use the values from the column with the short time format

Calulating sunrise and sunset in Python

I expected to find dozens of readily available implementations of sunrise and sunset calculations in Python on the web but this turned out to be a disappointment. Therefore I resolved to write my own straightforward implementation.

The sunrise equation is pretty easy to find on Wikipedia but actual implementations are not, certainly not in Python 3.x. (If you are willing to stick to Python 2.x there is of course the excellent PyEphem package, see the end of this article) Fortunately NOAA provides annotated equations in the form of an OpenOffice spreadsheet. This allows for simple re-engineering in Python and gives us a way to verify the results of a new implementation against those in the spreadsheet.

If you save the code below as sunrise.py then calculating the time of sunrise today would be straightforward as shown in this example:

 import datetime
 import sunrise
 s = sun(lat=49,long=3)
 print('sunrise at ',s.sunrise(when=datetime.datetime.now())

The sun class also provides a sunset() method and solarnoon() method. All three methods take a when parameter that should be a datetime.datetime object. If this object contains timezone information or daylight saving time information, this information is used when calculating the times of sunrise, sunset and the solar noon.

Note that if no when parameter is given, a default datetime is used that is initialized with a LocalTimezone object from the timezone module. I have not provided that module here but you can implement one simple enough by copying the example in Python's documentation or you can comment out the import statement below and always supply a when parameter.

from math import cos,sin,acos,asin,tan
from math import degrees as deg, radians as rad
from datetime import date,datetime,time

# this module is not provided here. See text.
from timezone import LocalTimezone

class sun:
 """ 
 Calculate sunrise and sunset based on equations from NOAA
 http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html

 typical use, calculating the sunrise at the present day:
 
 import datetime
 import sunrise
 s = sun(lat=49,long=3)
 print('sunrise at ',s.sunrise(when=datetime.datetime.now())
 """
 def __init__(self,lat=52.37,long=4.90): # default Amsterdam
  self.lat=lat
  self.long=long
  
 def sunrise(self,when=None):
  """
  return the time of sunrise as a datetime.time object
  when is a datetime.datetime object. If none is given
  a local time zone is assumed (including daylight saving
  if present)
  """
  if when is None : when = datetime.now(tz=LocalTimezone())
  self.__preptime(when)
  self.__calc()
  return sun.__timefromdecimalday(self.sunrise_t)
  
 def sunset(self,when=None):
  if when is None : when = datetime.now(tz=LocalTimezone())
  self.__preptime(when)
  self.__calc()
  return sun.__timefromdecimalday(self.sunset_t)
  
 def solarnoon(self,when=None):
  if when is None : when = datetime.now(tz=LocalTimezone())
  self.__preptime(when)
  self.__calc()
  return sun.__timefromdecimalday(self.solarnoon_t)
  
 @staticmethod
 def __timefromdecimalday(day):
  """
  returns a datetime.time object.
  
  day is a decimal day between 0.0 and 1.0, e.g. noon = 0.5
  """
  hours  = 24.0*day
  h      = int(hours)
  minutes= (hours-h)*60
  m      = int(minutes)
  seconds= (minutes-m)*60
  s      = int(seconds)
  return time(hour=h,minute=m,second=s)

 def __preptime(self,when):
  """
  Extract information in a suitable format from when, 
  a datetime.datetime object.
  """
  # datetime days are numbered in the Gregorian calendar
  # while the calculations from NOAA are distibuted as
  # OpenOffice spreadsheets with days numbered from
  # 1/1/1900. The difference are those numbers taken for 
  # 18/12/2010
  self.day = when.toordinal()-(734124-40529)
  t=when.time()
  self.time= (t.hour + t.minute/60.0 + t.second/3600.0)/24.0
  
  self.timezone=0
  offset=when.utcoffset()
  if not offset is None:
   self.timezone=offset.seconds/3600.0
  
 def __calc(self):
  """
  Perform the actual calculations for sunrise, sunset and
  a number of related quantities.
  
  The results are stored in the instance variables
  sunrise_t, sunset_t and solarnoon_t
  """
  timezone = self.timezone # in hours, east is positive
  longitude= self.long     # in decimal degrees, east is positive
  latitude = self.lat      # in decimal degrees, north is positive

  time  = self.time # percentage past midnight, i.e. noon  is 0.5
  day      = self.day     # daynumber 1=1/1/1900
 
  Jday     =day+2415018.5+time-timezone/24 # Julian day
  Jcent    =(Jday-2451545)/36525    # Julian century

  Manom    = 357.52911+Jcent*(35999.05029-0.0001537*Jcent)
  Mlong    = 280.46646+Jcent*(36000.76983+Jcent*0.0003032)%360
  Eccent   = 0.016708634-Jcent*(0.000042037+0.0001537*Jcent)
  Mobliq   = 23+(26+((21.448-Jcent*(46.815+Jcent*(0.00059-Jcent*0.001813))))/60)/60
  obliq    = Mobliq+0.00256*cos(rad(125.04-1934.136*Jcent))
  vary     = tan(rad(obliq/2))*tan(rad(obliq/2))
  Seqcent  = sin(rad(Manom))*(1.914602-Jcent*(0.004817+0.000014*Jcent))+sin(rad(2*Manom))*(0.019993-0.000101*Jcent)+sin(rad(3*Manom))*0.000289
  Struelong= Mlong+Seqcent
  Sapplong = Struelong-0.00569-0.00478*sin(rad(125.04-1934.136*Jcent))
  declination = deg(asin(sin(rad(obliq))*sin(rad(Sapplong))))
  
  eqtime   = 4*deg(vary*sin(2*rad(Mlong))-2*Eccent*sin(rad(Manom))+4*Eccent*vary*sin(rad(Manom))*cos(2*rad(Mlong))-0.5*vary*vary*sin(4*rad(Mlong))-1.25*Eccent*Eccent*sin(2*rad(Manom)))

  hourangle= deg(acos(cos(rad(90.833))/(cos(rad(latitude))*cos(rad(declination)))-tan(rad(latitude))*tan(rad(declination))))

  self.solarnoon_t=(720-4*longitude-eqtime+timezone*60)/1440
  self.sunrise_t  =self.solarnoon_t-hourangle*4/1440
  self.sunset_t   =self.solarnoon_t+hourangle*4/1440

if __name__ == "__main__":
 s=sun(lat=52.37,long=4.90)
 print(datetime.today())
 print(s.sunrise(),s.solarnoon(),s.sunset())

For people willing to stick to Python 2.x there is a simple and good alternative in the form of the PyEphem package. It can do a lot more than just calculating sunsise. An example is shown below.

import ephem
o=ephem.Observer()
o.lat='49'
o.long='3'
s=ephem.Sun()
s.compute()
print ephem.localtime(o.next_rising(s))