Saturday, 18 December 2010

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 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(

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

 typical use, calculating the sunrise at the present day:
 import datetime
 import sunrise
 s = sun(lat=49,long=3)
 print('sunrise at ',s.sunrise(
 def __init__(self,lat=52.37,long=4.90): # default Amsterdam
 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 =
  return sun.__timefromdecimalday(self.sunrise_t)
 def sunset(self,when=None):
  if when is None : when =
  return sun.__timefromdecimalday(self.sunset_t)
 def solarnoon(self,when=None):
  if when is None : when =
  return sun.__timefromdecimalday(self.solarnoon_t)
 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 = when.toordinal()-(734124-40529)
  self.time= (t.hour + t.minute/60.0 + t.second/3600.0)/24.0
  if not offset is None:
 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 =      # in decimal degrees, north is positive

  time  = self.time # percentage past midnight, i.e. noon  is 0.5
  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.sunrise_t  =self.solarnoon_t-hourangle*4/1440
  self.sunset_t   =self.solarnoon_t+hourangle*4/1440

if __name__ == "__main__":

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
print ephem.localtime(o.next_rising(s))


  1. I was all ready to write my own module in exactly the same fashion! Thanks so much for doing this :-)

  2. There is an error in the Western Hemisphere.
    line 80 should read:
    self.timezone=offset.seconds/3600.0 + (offset.days * 24)
    because time zones west of UTC have an offset value of -1 days plus a number of seconds.

  3. @VernonDCole You're right, thanks for the comment

  4. Why does this return UTC instead of local time? I'm getting sunset times of 24 hours, which chokes on. Why?

  5. @Anonymous: I don't understand your question. Can you give an example?

  6. Thanks for publishing this code. I took the liberty of reworking it to give just a day/night boolean return, and would appreciate your approval before releasing this revised code. Please see where we plan to use this for determining whether aircraft take off or land in day or night conditions.

    Many Thanks
    Dave Jesse

  7. Hi Michel, thanks for code and...anonymous, strangely I am just in process of doing the same thing :-)

  8. Thanks!
    Ephem has now an update for python 3.Xxx

  9. Hey Michel,

    cool script, thx...
    I used it in a pet project for my raspi.
    Is it ok for you to put it on GitHub?


  10. Hello,

    Thanks for the algorithm. When I try:

    sunrise.sun(49, -123).sunset(when=datetime.datetime(2013, 6, 1))

    I get:

    ValueError: hour must be in 0..23

    The error seems to occur for any longitude under -61. Any idea why this would occur?

    1. I'm also getting this error:

      s = sun(lat=-41.3,long=174.8)

    2. I get that error when the time zone does not match the longitude. Also, as a check I simply copied and pasted the formulas from the NOAA Excel spread sheet into Python (see below) and everything seems to work OK. Changes to the functions case (capitals to lower case) , removing the dollar sign from the absolute references (eg, $B$5 >> B5), and changing the MOD function from Excel format to Python format were the only changes required. The code then looked like this:


      # Calculate sunrise and sunset based on equations from NOAA

      # tested on IDLE

      from math import sin, asin, cos, acos, tan, atan2, radians, degrees

      def __timefromdecimalday(day):
      hours = 24.0*day
      h = int(hours)
      minutes= (hours-h)*60
      m = int(minutes)
      seconds= (minutes-m)*60
      s = int(seconds)
      return h,m,s

      time_zone = 10 # non DST
      lat = -38 # Melbourne, Australia
      lon = 145
      date = 41645 # number format from excel spreadsheet
      time = 0.1/24

      # setup parameters
      B3 = lat
      B4 = lon
      B5 = time_zone
      B7 = date
      # cut and pasted equations fron NOAA spreadsheet
      # may have to change case (capitals to lower case, etc) of functions
      D2 = B7
      E2 = 0.1/24 # day is a decimal fraction of 24 hours between 0.0 and 1.0 (e.g. noon = 0.5)
      F2 = D2 + 2415018.5 + E2 - B5 / 24
      G2 = (F2 - 2451545) / 36525
      I2 = (280.46646+G2*(36000.76983 + G2*0.0003032)) % 360
      J2 = 357.52911+G2*(35999.05029 - 0.0001537*G2)
      K2 = 0.016708634-G2*(0.000042037+0.0000001267*G2)
      L2 = sin(radians(J2))*(1.914602-G2*(0.004817+0.000014*G2))+sin(radians(2*J2))*(0.019993-0.000101*G2)+sin(radians(3*J2))*0.000289
      M2 = I2 + L2
      P2 = M2-0.00569-0.00478*sin(radians(125.04-1934.136*G2))
      Q2 = 23+(26+((21.448-G2*(46.815+G2*(0.00059-G2*0.001813))))/60)/60
      R2 = Q2+0.00256*cos(radians(125.04-1934.136*G2))
      T2 = degrees(asin(sin(radians(R2))*sin(radians(P2))))
      U2 = tan(radians(R2/2))*tan(radians(R2/2))
      V2 = 4*degrees(U2*sin(2*radians(I2))-2*K2*sin(radians(J2))+4*K2*U2*sin(radians(J2))*cos(2*radians(I2))-0.5*U2*U2*sin(4*radians(I2))-1.25*K2*K2*sin(2*radians(J2)))
      W2 = degrees(acos(cos(radians(90.833))/(cos(radians(B3))*cos(radians(T2)))-tan(radians(B3))*tan(radians(T2))))
      X2 = (720-4*B4-V2+B5*60)/1440
      Y2 = X2-W2*4/1440 # sunrise
      Z2 = X2+W2*4/1440 # sunset

      hrs,mins,secs = __timefromdecimalday(Y2)
      print "sunrise {:02}:{:02}:{:02}".format(hrs,mins,secs)

      hrs,mins,secs = __timefromdecimalday(Z2)
      print "sunset {:02}:{:02}:{:02}".format(hrs,mins,secs)

  11. Sorry, i´m a absolute beginner in Python and programming. The works fine, but when i try the typical use, i get an syntax error in the line "print('sunrise at ',s.sunrise( "
    if i comment this line out, i get the error " name sun is not defined"...
    whats wrong?

    1. I ran into the same issue using Python 3 on a Win7 box.
      import datetime
      import sunrise

      # Latitude N is positive, Latitude S is negative
      # Longitude E is positive, Longitude W is negative

      s = sunrise.sun(lat=25,long=-80)

  12. here is my solution:

    import calendar
    from math import cos, sin, acos as arccos, asin as arcsin, tan as tg, degrees, radians

    def mod(a,b):
    return a % b

    def isLeapYear(year):
    return (year % 4 == 0 and year % 100 != 0) or year % 400 == 0

    def getDayNumber(year, month, day):
    cnt = 0
    for t in range(1900,year):
    if isLeapYear(t):
    cnt += 366
    cnt += 365
    for t in range(1,month):
    cnt += calendar.monthrange(2014, 2)[1]
    return cnt + day + 1

    def getHHMMSS(h):
    hh = int(h)
    mm = (h - hh) * 60
    ss = int(0.5 + (mm - int(mm)) * 60)
    return "{0:2d}:{1:02d}:{2:02d}" . format(hh, int(mm), ss)

    # based on:
    def getSunriseAndSunset(lat, lon, dst, year, month, day):
    localtime = 12.00
    b2 = lat
    b3 = lon
    b4 = dst
    b5 = localtime / 24
    b6 = year
    d30 = getDayNumber(year, month, day)
    e30 = b5
    f30 = d30 + 2415018.5 + e30 - b4 / 24
    g30 = (f30 - 2451545) / 36525
    q30 = 23 + (26 + ((21.448 - g30 * (46.815 + g30 * (0.00059 - g30 * 0.001813)))) / 60) / 60
    r30 = q30 + 0.00256 * cos(radians(125.04 - 1934.136 * g30))
    j30 = 357.52911 + g30 * (35999.05029 - 0.0001537 * g30)
    k30 = 0.016708634 - g30 * (0.000042037 + 0.0000001267 * g30)
    l30 = sin(radians(j30)) * (1.914602 - g30 * (0.004817 + 0.000014 * g30)) + sin(radians(2 * j30)) * (0.019993 - 0.000101 * g30) + sin(radians(3 * j30)) * 0.000289
    i30 = mod(280.46646 + g30 * (36000.76983 + g30 * 0.0003032), 360)
    m30 = i30 + l30
    p30 = m30 - 0.00569 - 0.00478 * sin(radians(125.04 - 1934.136 * g30))
    t30 = degrees(arcsin(sin(radians(r30)) * sin(radians(p30))))
    u30 = tg(radians(r30 / 2)) * tg(radians(r30 / 2))
    v30 = 4 * degrees(u30 * sin(2 * radians(i30)) - 2 * k30 * sin(radians(j30)) + 4 * k30 * u30 * sin(radians(j30)) * cos(2 * radians(i30)) - 0.5 * u30 * u30 * sin(4 * radians(i30)) - 1.25 * k30 * k30 * sin(2 * radians(j30)))
    w30 = degrees(arccos(cos(radians(90.833)) / (cos(radians(b2)) * cos(radians(t30))) - tg(radians(b2)) * tg(radians(t30))))
    x30 = (720 - 4 * b3 - v30 + b4 * 60) / 1440
    x30 = (720 - 4 * b3 - v30 + b4 * 60) / 1440
    y30 = (x30 * 1440 - w30 * 4) / 1440
    z30 = (x30 * 1440 + w30 * 4) / 1440
    sunrise = y30 * 24
    sunset = z30 * 24
    return (sunrise, sunset)

    # latitude (+N -S):
    lat = 50.0877777777777
    # longitude (+E -W):
    lon = 14.4205555555555
    # (+E -W)
    dst = 1
    (sunrise, sunset) = getSunriseAndSunset(lat, lon, dst, 2014, 1, 29)
    print("sunrise=", getHHMMSS(sunrise))
    print("sunset =", getHHMMSS(sunset))

    1. please change line
      cnt += calendar.monthrange(2014, 2)[1]
      to this
      cnt += calendar.monthrange(year, t)[1]

    2. As far as I see Your x30-line is double . . .

  13. Thanks a lot for sharing this module. I will be going to use it in my selfmade home automation system to control the electric shutters.

  14. Thanks for sharing. This made my day easier.

  15. Unfortunately you have not run into Python Linters yet. These are not annoying rules that critic our code but ways to help us make it more readable. From my experience you have some functions with way too many lines of code. Python likes at least three letter variables and less than 25 lines of code in a function. But that's these days not Py2.7 2010

  16. Replies
    1. same here.
      timezone or timezones? But there is no LocalTimezone object in the timezones

    2. You have to provide an object yourself for that.
      If you use Python 3.7 I have made modified a version that does that:

  17. Hello, there is a mistake in the Eccent formula. In the python code, the coef for J^2 is 0.0001537 (same as Manom).
    The correct coef is 0.0000001267 (XLS file from NOAA, column K).

  18. Thanks, this is amazing!
    I'd like to use it to create a CLI tool, of course mentioning the source of the original code. Is it OK?

  19. This professional hacker is absolutely reliable and I strongly recommend him for any type of hack you require. I know this because I have hired him severally for various hacks and he has never disappointed me nor any of my friends who have hired him too, he can help you with any of the following hacks:

    -Phone hacks (remotely)
    -Credit repair
    -Bitcoin recovery (any cryptocurrency)
    -Make money from home (USA only)
    -Social media hacks
    -Website hacks
    -Erase criminal records (USA & Canada only)
    -Grade change
    -funds recovery

    Email: onlineghosthacker247@ gmail .com

  20. Actually I read it yesterday I looked at most of your posts but I had some ideas about it . This article is probably where I got the most useful information for my research and today I wanted to read it again because it is so well written.
    Data Science Course in Bangalore

  21. Top quality blog with unique content and information shared was valuable looking forward for next updated thank you.
    Ethical Hacking Course in Bangalore