I want to calculate the sunrise and sunset times for a full year for a given location taking topography into account. Perhaps sunrise/sunset are not the right terms, but what I want is the time at which the sun rises above the horizon, and the time at which it drops below the horizon, taking any hills into account.
I am using QGIS or GRASS. I can use r.horizon to generate the horizon angle from a specific point, but I'm not sure how to get from there to the sunrise/sunset times.
I installed the ephem python package in my GNU/Debian Linux System and, I could use it at the Python Console of QGIS. I created an observer for a point near of Utah Lake (USA) and it was calculated the raising and setting sun for '2010/6/21' day. The script was:
import ephem
#defining an observer
obs = ephem.Observer()
#defining position
long = '-112.092807'
lat = '40.135114'
obs.long = ephem.degrees(long)
obs.lat = ephem.degrees(lat)
print "long = ", obs.long, "lat = ", obs.lat
#defining date
date = '2010/6/21'
obs.date = ephem.Date(date)
#defining an astronomic object; Sun in this case
sun = ephem.Sun(obs)
r1 = obs.next_rising(sun)
s1 = obs.next_setting(sun)
print "rising sun (UTC time): ", r1
print "setting sun (UTC time): ", s1
r1_lt = ephem.Date(r1 - 6 * ephem.hour) #local time
(y, mn, d, h, min, s) = r1_lt.tuple()
print "rising sun: (local time): {:.2f}".format( h + min/60. + s/3600. )
s1_lt = ephem.Date(s1 - 6 * ephem.hour) #local time
(y, mn, d, h, min, s) = s1_lt.tuple()
print "setting sun (local time): {:.2f}".format( h + min/60. + s/3600. )
After running at the Python Console of QGIS, the result was:
>>>execfile(u'/home/zeito/pyqgis_scripts/ephem.py'.encode('UTF-8'))
long = -112:05:34.1 lat = 40:08:06.4
rising sun (UTC time): 2010/6/21 11:58:58
setting sun (UTC time): 2010/6/21 03:01:14
rising sun: (local time): 5.98
setting sun (local time): 21.02
This is the answer.
Editing note:
Defining a new horizon (e.g. 5 degrees):
.
.
.
obs.horizon = '5'
sun = ephem.Sun(obs)
r1 = obs.next_rising(sun)
s1 = obs.next_setting(sun)
print "rising sun (UTC time): ", r1
print "setting sun (UTC time): ", s1
r1_lt = ephem.Date(r1 - 6 * ephem.hour) #local time
(y, mn, d, h, min, s) = r1_lt.tuple()
print "rising sun: (local time): {:.2f}".format( h + min/60. + s/3600. )
s1_lt = ephem.Date(s1 - 6 * ephem.hour) #local time
(y, mn, d, h, min, s) = s1_lt.tuple()
print "setting sun (local time): {:.2f}".format( h + min/60. + s/3600. )
the result is:
>>>execfile(u'/home/zeito/pyqgis_scripts/ephem.py'.encode('UTF-8'))
long = -112:05:34.1 lat = 40:08:06.4
rising sun (UTC time): 2010/6/21 12:31:48
setting sun (UTC time): 2010/6/21 02:28:24
rising sun: (local time): 6.53
setting sun (local time): 20.47
You could use a library like PyEphem to find, for a given location, elevation, day, time(s) (and planet ;) the azimut and altitude angles of the sun.
You can use a little bit of python to do this but you need first to generate a dictionary or table, from your observers position, of the angle to the real world horizon seems to be at. This would need to be at say 1 degree intervals with 0 = true north and for each degree giving the angle of the horizon. This can come from your r.horizon.
You could then, use pyephem to a) create an observer at the latitude, longitude and altitude of your observer and for each minute of a given day use it to calculate the suns apparent positions, these would be given as four values, ra, dec, alt and az.
Then for each time value you can compare the alt with the value from the elevation table at that azimuth and when the alt is greater than the value in the table you will be able to see the sun.
It is even possible that you can find some points where, on a given day, the sun rises in a valley and then passes out of site behind a mountain then re-appears over or the other side of the mountain.
