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:

view plainprint?
  1. import datetime  
  2. import sunrise  
  3. s = sun(lat=49,long=3)  
  4. 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.

view plainprint?
  1. from math import cos,sin,acos,asin,tan  
  2. from math import degrees as deg, radians as rad  
  3. from datetime import date,datetime,time  
  4.   
  5. # this module is not provided here. See text.  
  6. from timezone import LocalTimezone  
  7.   
  8. class sun:  
  9.  """  
  10.  Calculate sunrise and sunset based on equations from NOAA 
  11.  http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html 
  12.  
  13.  typical use, calculating the sunrise at the present day: 
  14.   
  15.  import datetime 
  16.  import sunrise 
  17.  s = sun(lat=49,long=3) 
  18.  print('sunrise at ',s.sunrise(when=datetime.datetime.now()) 
  19.  """  
  20.  def __init__(self,lat=52.37,long=4.90): # default Amsterdam  
  21.   self.lat=lat  
  22.   self.long=long  
  23.     
  24.  def sunrise(self,when=None):  
  25.   """ 
  26.   return the time of sunrise as a datetime.time object 
  27.   when is a datetime.datetime object. If none is given 
  28.   a local time zone is assumed (including daylight saving 
  29.   if present) 
  30.   """  
  31.   if when is None : when = datetime.now(tz=LocalTimezone())  
  32.   self.__preptime(when)  
  33.   self.__calc()  
  34.   return sun.__timefromdecimalday(self.sunrise_t)  
  35.     
  36.  def sunset(self,when=None):  
  37.   if when is None : when = datetime.now(tz=LocalTimezone())  
  38.   self.__preptime(when)  
  39.   self.__calc()  
  40.   return sun.__timefromdecimalday(self.sunset_t)  
  41.     
  42.  def solarnoon(self,when=None):  
  43.   if when is None : when = datetime.now(tz=LocalTimezone())  
  44.   self.__preptime(when)  
  45.   self.__calc()  
  46.   return sun.__timefromdecimalday(self.solarnoon_t)  
  47.    
  48.  @staticmethod  
  49.  def __timefromdecimalday(day):  
  50.   """ 
  51.   returns a datetime.time object. 
  52.    
  53.   day is a decimal day between 0.0 and 1.0, e.g. noon = 0.5 
  54.   """  
  55.   hours  = 24.0*day  
  56.   h      = int(hours)  
  57.   minutes= (hours-h)*60  
  58.   m      = int(minutes)  
  59.   seconds= (minutes-m)*60  
  60.   s      = int(seconds)  
  61.   return time(hour=h,minute=m,second=s)  
  62.   
  63.  def __preptime(self,when):  
  64.   """ 
  65.   Extract information in a suitable format from when,  
  66.   a datetime.datetime object. 
  67.   """  
  68.   # datetime days are numbered in the Gregorian calendar  
  69.   # while the calculations from NOAA are distibuted as  
  70.   # OpenOffice spreadsheets with days numbered from  
  71.   # 1/1/1900. The difference are those numbers taken for   
  72.   # 18/12/2010  
  73.   self.day = when.toordinal()-(734124-40529)  
  74.   t=when.time()  
  75.   self.time= (t.hour + t.minute/60.0 + t.second/3600.0)/24.0  
  76.     
  77.   self.timezone=0  
  78.   offset=when.utcoffset()  
  79.   if not offset is None:  
  80.    self.timezone=offset.seconds/3600.0  
  81.     
  82.  def __calc(self):  
  83.   """ 
  84.   Perform the actual calculations for sunrise, sunset and 
  85.   a number of related quantities. 
  86.    
  87.   The results are stored in the instance variables 
  88.   sunrise_t, sunset_t and solarnoon_t 
  89.   """  
  90.   timezone = self.timezone # in hours, east is positive  
  91.   longitude= self.long     # in decimal degrees, east is positive  
  92.   latitude = self.lat      # in decimal degrees, north is positive  
  93.   
  94.   time  = self.time # percentage past midnight, i.e. noon  is 0.5  
  95.   day      = self.day     # daynumber 1=1/1/1900  
  96.    
  97.   Jday     =day+2415018.5+time-timezone/24 # Julian day  
  98.   Jcent    =(Jday-2451545)/36525    # Julian century  
  99.   
  100.   Manom    = 357.52911+Jcent*(35999.05029-0.0001537*Jcent)  
  101.   Mlong    = 280.46646+Jcent*(36000.76983+Jcent*0.0003032)%360  
  102.   Eccent   = 0.016708634-Jcent*(0.000042037+0.0001537*Jcent)  
  103.   Mobliq   = 23+(26+((21.448-Jcent*(46.815+Jcent*(0.00059-Jcent*0.001813))))/60)/60  
  104.   obliq    = Mobliq+0.00256*cos(rad(125.04-1934.136*Jcent))  
  105.   vary     = tan(rad(obliq/2))*tan(rad(obliq/2))  
  106.   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  
  107.   Struelong= Mlong+Seqcent  
  108.   Sapplong = Struelong-0.00569-0.00478*sin(rad(125.04-1934.136*Jcent))  
  109.   declination = deg(asin(sin(rad(obliq))*sin(rad(Sapplong))))  
  110.     
  111.   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)))  
  112.   
  113.   hourangle= deg(acos(cos(rad(90.833))/(cos(rad(latitude))*cos(rad(declination)))-tan(rad(latitude))*tan(rad(declination))))  
  114.   
  115.   self.solarnoon_t=(720-4*longitude-eqtime+timezone*60)/1440  
  116.   self.sunrise_t  =self.solarnoon_t-hourangle*4/1440  
  117.   self.sunset_t   =self.solarnoon_t+hourangle*4/1440  
  118.   
  119. if __name__ == "__main__":  
  120.  s=sun(lat=52.37,long=4.90)  
  121.  print(datetime.today())  
  122.  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.

view plainprint?
  1. import ephem  
  2. o=ephem.Observer()  
  3. o.lat='49'  
  4. o.long='3'  
  5. s=ephem.Sun()  
  6. s.compute()  
  7. print ephem.localtime(o.next_rising(s))  
Labels: , , ,