11

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.

  • Sounds like a reverse viewshed... at what angle can the observer (sun) see the location then use that as an offset to the standard sunrise/set times. That's a very good question. – Michael Stimson Apr 26 '16 at 0:50
2

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
  • Thank you for this, it seems like a piece of the puzzle. It seems like I could use PyEphem to calculate the azimuth and altitude angles of the sun at a particular time? How do I then use that to find the time when the sun rises above/sets below the horizon? (Accounting for topography, I'm assuming that the sunset/sunrise times that the ephem package is outputting are assuming a perfectly smooth earth.) – Stu Apr 26 '16 at 16:11
  • 2
    Perhaps some combination of r.horizon and PyEphem would work? Any thoughts on how to use the two together? GRASS does have a module r.sun, which I can use to calculate the number of hours of direct sunlight at a given location on a given day, accounting for topography. That seems to almost do what I want, but it doesn’t seem to output the actual time the direct sunlight starts/stops. – Stu Apr 26 '16 at 16:27
  • Promising and useful, but it does not take into account topography. The visible horizon from the inside of a valley is not the same as the true horizon (unobscured). – alphabetasoup Apr 26 '16 at 20:02
  • @RichardLaw You can set the topography with the 'horizon' method. – xunilk Apr 27 '16 at 6:04
  • 1
    It seems like the horizon method is just setting a single horizon angle for the entire horizon? If we take topography into account, the horizon is not perfectly flat. So, the angle to the horizon will be different for each different azimuth that is calculated. The r.horizon module will output the angle to the horizon for each different azimuth. It won't be a single number (eg. 5 degrees), but rather a data set of a series of elevations for all azimuths (5 degrees elevation at 180 degrees azimuth, 6 degrees elevation at 185 degrees azimuth, 7 degrees elevation at 190 degrees azimuth, etc.) – Stu Apr 27 '16 at 19:54
1

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.

  • Can you suggest an implementation of this? The original question makes no mention of python implicitly but the O.P. might be able to pick it up if the example is simple enough. – Michael Stimson Apr 26 '16 at 4:43
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.