1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
|
import calendar
import math
import datetime
from dateutil import tz
class SunTimeException(Exception):
def __init__(self, message):
super(SunTimeException, self).__init__(message)
class Sun:
"""
Approximated calculation of sunrise and sunset datetimes. Adapted from:
https://stackoverflow.com/questions/19615350/calculate-sunrise-and-sunset-times-for-a-given-gps-coordinate-within-postgresql
"""
def __init__(self, lat, lon):
self._lat = lat
self._lon = lon
def get_sunrise_time(self, date=None):
"""
Calculate the sunrise time for given date.
:param lat: Latitude
:param lon: Longitude
:param date: Reference date. Today if not provided.
:return: UTC sunrise datetime
:raises: SunTimeException when there is no sunrise and sunset on given location and date
"""
date = datetime.date.today() if date is None else date
sr = self._calc_sun_time(date, True)
if sr is None:
raise SunTimeException('The sun never rises on this location (on the specified date)')
else:
return sr
def get_local_sunrise_time(self, date=None, local_time_zone=tz.tzlocal()):
"""
Get sunrise time for local or custom time zone.
:param date: Reference date. Today if not provided.
:param local_time_zone: Local or custom time zone.
:return: Local time zone sunrise datetime
"""
date = datetime.date.today() if date is None else date
sr = self._calc_sun_time(date, True)
if sr is None:
raise SunTimeException('The sun never rises on this location (on the specified date)')
else:
return sr.astimezone(local_time_zone)
def get_sunset_time(self, date=None):
"""
Calculate the sunset time for given date.
:param lat: Latitude
:param lon: Longitude
:param date: Reference date. Today if not provided.
:return: UTC sunset datetime
:raises: SunTimeException when there is no sunrise and sunset on given location and date.
"""
date = datetime.date.today() if date is None else date
ss = self._calc_sun_time(date, False)
if ss is None:
raise SunTimeException('The sun never sets on this location (on the specified date)')
else:
return ss
def get_local_sunset_time(self, date=None, local_time_zone=tz.tzlocal()):
"""
Get sunset time for local or custom time zone.
:param date: Reference date
:param local_time_zone: Local or custom time zone.
:return: Local time zone sunset datetime
"""
date = datetime.date.today() if date is None else date
ss = self._calc_sun_time(date, False)
if ss is None:
raise SunTimeException('The sun never sets on this location (on the specified date)')
else:
return ss.astimezone(local_time_zone)
def _calc_sun_time(self, date, isRiseTime=True, zenith=90.8):
"""
Calculate sunrise or sunset date.
:param date: Reference date
:param isRiseTime: True if you want to calculate sunrise time.
:param zenith: Sun reference zenith
:return: UTC sunset or sunrise datetime
:raises: SunTimeException when there is no sunrise and sunset on given location and date
"""
# isRiseTime == False, returns sunsetTime
day = date.day
month = date.month
year = date.year
TO_RAD = math.pi/180.0
# 1. first calculate the day of the year
N1 = math.floor(275 * month / 9)
N2 = math.floor((month + 9) / 12)
N3 = (1 + math.floor((year - 4 * math.floor(year / 4) + 2) / 3))
N = N1 - (N2 * N3) + day - 30
# 2. convert the longitude to hour value and calculate an approximate time
lngHour = self._lon / 15
if isRiseTime:
t = N + ((6 - lngHour) / 24)
else: #sunset
t = N + ((18 - lngHour) / 24)
# 3. calculate the Sun's mean anomaly
M = (0.9856 * t) - 3.289
# 4. calculate the Sun's true longitude
L = M + (1.916 * math.sin(TO_RAD*M)) + (0.020 * math.sin(TO_RAD * 2 * M)) + 282.634
L = self._force_range(L, 360 ) #NOTE: L adjusted into the range [0,360)
# 5a. calculate the Sun's right ascension
RA = (1/TO_RAD) * math.atan(0.91764 * math.tan(TO_RAD*L))
RA = self._force_range(RA, 360 ) #NOTE: RA adjusted into the range [0,360)
# 5b. right ascension value needs to be in the same quadrant as L
Lquadrant = (math.floor( L/90)) * 90
RAquadrant = (math.floor(RA/90)) * 90
RA = RA + (Lquadrant - RAquadrant)
# 5c. right ascension value needs to be converted into hours
RA = RA / 15
# 6. calculate the Sun's declination
sinDec = 0.39782 * math.sin(TO_RAD*L)
cosDec = math.cos(math.asin(sinDec))
# 7a. calculate the Sun's local hour angle
cosH = (math.cos(TO_RAD*zenith) - (sinDec * math.sin(TO_RAD*self._lat))) / (cosDec * math.cos(TO_RAD*self._lat))
if cosH > 1:
return None # The sun never rises on this location (on the specified date)
if cosH < -1:
return None # The sun never sets on this location (on the specified date)
# 7b. finish calculating H and convert into hours
if isRiseTime:
H = 360 - (1/TO_RAD) * math.acos(cosH)
else: #setting
H = (1/TO_RAD) * math.acos(cosH)
H = H / 15
#8. calculate local mean time of rising/setting
T = H + RA - (0.06571 * t) - 6.622
#9. adjust back to UTC
UT = T - lngHour
UT = self._force_range(UT, 24) # UTC time in decimal format (e.g. 23.23)
#10. Return
hr = self._force_range(int(UT), 24)
min = round((UT - int(UT))*60, 0)
if min == 60:
hr += 1
min = 0
#10. check corner case https://github.com/SatAgro/suntime/issues/1
if hr == 24:
hr = 0
day += 1
if day > calendar.monthrange(year, month)[1]:
day = 1
month += 1
if month > 12:
month = 1
year += 1
return datetime.datetime(year, month, day, hr, int(min), tzinfo=tz.tzutc())
@staticmethod
def _force_range(v, max):
# force v to be >= 0 and < max
if v < 0:
return v + max
elif v >= max:
return v - max
return v
if __name__ == '__main__':
sun = Sun(85.0, 21.00)
try:
print(sun.get_local_sunrise_time())
print(sun.get_local_sunset_time())
# On a special date in UTC
abd = datetime.date(2014, 1, 3)
abd_sr = sun.get_local_sunrise_time(abd)
abd_ss = sun.get_local_sunset_time(abd)
print(abd_sr)
print(abd_ss)
except SunTimeException as e:
print("Error: {0}".format(e))
|
