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
sunclass also provides a
solarnoon()method. All three methods take a
whenparameter that should be a
datetime.datetimeobject. 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
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))