Observer.sun_set_time not working properly

[vlizanae]

I'm trying to use Observer to determine the time of the sunset for a given day, the problem is that here on the site we have an official sunset time that does not match the one that astroplan gives. After some observation the official time seems to be the correct, the times differ for about 10 minutes.

What worries me the most is that while debugging I realized that every time the code was executed the sunset time returned was different.

Minimal example:

from astroplan import Observer
from astropy.time import Time
from astropy.units import deg, m
from astropy.coordinates import EarthLocation
from datetime import datetime

obs = Observer(latitude=-24.622997508*deg, longitude=-70.40249839*deg, elevation=2635*m, name='pl', timezone='UTC')
print(obs.sun_set_time(Time(datetime.now()), which='next').to_datetime())

[bmorris3]

Hi @vlizanae! Thanks for this issue, let's see if we can make some sense of it, in two parts:

Matching the official time

One source of difference is likely due to the way that you initialized the Observer object. Observer takes an optional keyword called pressure which sets the atmospheric pressure where your observations are occurring. By default the apparent position sky transformations are done without an atmosphere unless pressure is not None. Try setting

from astropy.units import bar
obs = Observer(latitude=-24.622997508*deg, longitude=-70.40249839*deg, 
               elevation=2635*m, name='pl', timezone='UTC', pressure=1*bar)
print(obs.sun_set_time(Time(datetime.now()), which='next').to_datetime())

and see how that compares to your official sunset time.

---

Changing sunset time

The change in the sunset time when making repeated calls to the function arises from supplying the function with datetime.now(), which is constantly changing, in combination with the approximation we make to calculate target rise/set times. In summary, we create a grid of times over the next/previous 24 hours and look for horizon crossings of the target. Precisely where the grid points fall depends, and therefore the precision on the rise time (on the order of seconds), depends on the time you initialize on.

If instead you did this:

now = datetime.now()
print(obs.sun_set_time(Time(now), which='next').to_datetime())
print(obs.sun_set_time(Time(now), which='next').to_datetime())

you should see that the same sunset time should be returned as long as now is a constant.

[vlizanae]

Thanks for the quick response,

No problem with the approximation, I just wasn't expecting the given time to affect the time of the next sunset.

I tried with 1 bar and also by querying the value of the pressure from the same system that reports the sunset time but the result is more or less the same, the pressure needed to match the system's time is completely off the charts.

[bmorris3]

Interesting – another thing which might affect our calculations is that astroplan is computing the time that the center of the sun rises above the horizon – do you know if your system reports the rise of the solar limb?

[vlizanae]

It doesn't have a lot of information on the sun actually, it is mostly for other variables.

[bmorris3]

May I ask what this other software this is? Perhaps we can find out more details about it.

Our sunrise/set functions are validated against other astronomical packages and the tests are still passing, so it'd be interesting to know what this other software is doing under the hood.

[bmorris3]

Hi @vlizanae – is this still a problem for you, or can we close this issue? Thanks!

[tepickering]

i'm seeing this as well. astroplan gives sunset times that are consistently ~3-4 min earlier than what's published in our almanac, http://www.mmto.org/sites/default/files/almanac_2019.pdf. the almanac times have been verified by telescope operators for many years and are consistent with what i get with various online calculators like the one at NOAA, https://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html.

[bmorris3]

Thanks @tepickering. Is it possible this could be related to: https://github.com/astropy/astroplan/issues/409#issuecomment-500952670?

[tepickering]

i don't think so. i tried different horizon levels between +/- 1 deg and astroplan still predicts significantly earlier sunset time. i also noticed that the test against pyephem is done at the equator. would be worth comparing at 1 or 2 different sites at different latitudes. it's too bad the USNO site is down because that's kind of the ground truth for these kinds of calculations.

[tepickering]

ok, i think i've solved the riddle, at least to my satisfaction (fwiw). turns out the USNO definition for sunset/sunrise is when the center of the sun is 0.8333 degrees below the horizon (see https://rhodesmill.org/skyfield/almanac.html#sunrise-and-sunset while the USNO page is offline). this is an approximation that takes into account the sun's average radius and an average amount of atmospheric refraction. so to compare astroplan's times with other sources you need to set pressure=0*u.mbar in Observer and horizon=-0.8333*u.deg in Observer.sun_set_time() and Observer.sun_rise_time(). if i do this, then the astroplan rise/set times agree with skyfield's times to within a second and are in agreement with our almanac (which rounds to the nearest minute).

[bmorris3]

@tepickering Thanks so much for your superlative sunset sleuthing! Would you be interested in collaborating to add some tests to astroplan that validate the astroplan rise/set functions in these cases? I'll write up a PR that explains this caveat in the narrative docs, and validates against the MMTO almanac to the nearest minute, your comments there would be appreciated :+1:.