Wrong calculation of obscuring feature calculation

[Nechoj]

I think there is a bug in method adjust_to_obscuring_feature in module sun.py in line 249:

If elevation[0] is the height difference and elevation[1] is the distance to obscuring feature, then the calculation should use asin instead of acos:

asin(fabs(elevation[0]) / sqrt(pow(elevation[0], 2) + pow(elevation[1], 2)))

[pnbruckner]

I think it's more than that. The initial test of elevation[0] == 0.0 should be elevation[1] == 0.0, and the resulting sign is wrong. This value is added to the "zenith", so should be negative for positive heights and vice versa.

I'm stuck on version 2.2, and my workaround is to negate distance and reverse the elements of the tuple, so pass: (-distance, height).

[pnbruckner]

Actually, that only works for positive heights. So, to work for both positive and negative (relative) heights, the workaround should be:

observer_elevation = (-copysign(1, height) * distance, fabs(height))

Oh, and apparently, it's possible for floats to be -0.0 in addition to 0.0 where -0.0 != 0.0. So, that initial test should probably also be elevation[1] in (0.0, -0.0), or maybe just do away with it, since 1) the result will be the same, and 2) comparing floats for equality is not always a good thing to depend on.

[pnbruckner]

Here is what I'd recommend:

def adjust_to_obscuring_feature(elevation: Tuple[float, float]) -> float:
    """Calculate the number of degrees to adjust for an obscuring feature"""
    height, distance = elevation
    assert not distance < 0
    return -copysign(1, height) * degrees(
        acos(distance / sqrt(distance * distance + height * height))
    )

[Nechoj]

Please note that the parameter elevation defines the position of the observer relative to the obstacle and not the geometry of the obstacle (see also https://sffjunkie.github.io/astral/index.html#elevation-of-the-sun). E.g.: elevation=(300, 5000) means that the observer is elevated by 300m to the object/horizon at 5000m distance, where the object defines the horizon for surise and sunset. Likewise, elevation=(-300, 5000) means that the observer is 300m below the object/horizon at 5000m distance.

Test cases for the correct implementation therefore include:
elevation=(0, 5000) => times for sunrise and sunset remain unchanged compared to calculation without elevation
elevation=(10, 5000) => time for sunrise slightly earlier and time for sunset slightly later
elevation=(500, 5000) => time for sunrise significantly earlier and time for sunset significantly later
elevation=(-10, 5000) => time for sunrise slightly later and time for sunset slightly earlier
elevation=(-500, 5000) => time for sunrise significantly later and time for sunset significantly earlier

Given that, the correct implementation is achieved by replacing acos with asin as suggested in the first comment:

def adjust_to_obscuring_feature(elevation: Tuple[float, float]) -> float:
    """ Calculate the number of degrees to adjust for an obscuring feature"""
    if elevation[0] == 0.0:
        return 0.0

    sign = -1 if elevation[0] < 0.0 else 1
    return sign * math.degrees(
        math.asin(math.fabs(elevation[0]) / math.sqrt(pow(elevation[0], 2) + pow(elevation[1], 2)))
    )

(the other parts of the function are taken from https://github.com/sffjunkie/astral/blob/master/src/astral/sun.py without modification. Especially, elevation[0]=0 is a valid input and should return 0.0)

[pnbruckner]

Your implementation produces the same values as what I suggested, but still with an incorrect sign.

>>> for e in ((0, 5000), (10, 5000), (500, 5000), (-10, 5000), (-500, 5000)):
...     print(adjust_to_obscuring_feature_mine(e), adjust_to_obscuring_feature(e))
... 
-0.0 0.0
-0.11459140623948247 0.11459140623778596
-5.710593137499696 5.710593137499642
0.11459140623948247 -0.11459140623778596
5.710593137499696 -5.710593137499642

E.g., you say:

elevation=(10, 5000) => time for sunrise slightly earlier and time for sunset slightly later

If the obstruction is 10 meters higher than the observer, then the sunrise should be slightly later, and the sunset should be slightly earlier.

Also, "Especially, elevation[0]=0 is a valid input and should return 0.0)", which is true, but you see it still does even without that:

>>> def adjust_to_obscuring_feature(elevation: Tuple[float, float]) -> float:
...     """ Calculate the number of degrees to adjust for an obscuring feature"""
...     sign = -1 if elevation[0] < 0.0 else 1
...     return sign * math.degrees(
...         math.asin(math.fabs(elevation[0]) / math.sqrt(pow(elevation[0], 2) + pow(elevation[1], 2)))
...     )
... 
>>> adjust_to_obscuring_feature((0, 5000))
0.0

So, bottom line, either way works, as long as yours is changed to negate the return value.

[pnbruckner]

Please note that the parameter elevation defines the position of the observer relative to the obstacle and not the geometry of the obstacle (see also https://sffjunkie.github.io/astral/index.html#elevation-of-the-sun). E.g.: elevation=(300, 5000) means that the observer is elevated by 300m to the object/horizon at 5000m distance, where the object defines the horizon for surise and sunset.

I think you quoted the wrong part of the docs. You should look here: https://sffjunkie.github.io/astral/index.html#times-of-the-sun

Especially where it says, "For the second case pass a tuple of 2 floats. The first being the vertical distance to the top of the feature and the second the horizontal distance to the feature." This clearly means the opposite of what you explained. Also, although poorly documented in the code, there are places where it says the same thing.

I think it's reasonable to state the documentation and implementation are currently in need of change with respect to this issue.

Which makes me wonder. Is this package still being maintained? There hasn't been a release or commit in over a year.

[pnbruckner]

FWIW:

>>> timeit.timeit("adjust_to_obscuring_feature_mine((500, 5000))", globals=globals())
0.3826604629866779
>>> timeit.timeit("adjust_to_obscuring_feature((500, 5000))", globals=globals())
0.47058659396134317

[Nechoj]

You did not test the module as a whole calculating sunset and sunrise. The return value of adjust_to_obscuring_feature must meet the requirements of the other methods in the module. Example:

import astral
import datetime as dt

def adjust_to_obscuring_feature(elevation: Tuple[float, float]) -> float:
    """ Calculate the number of degrees to adjust for an obscuring feature
    correcting error in astral.sun: https://github.com/sffjunkie/astral/blob/master/src/astral/sun.py line 249
    use asin instead of acos
    """

    if elevation[0] == 0.0:
        return 0.0

    sign = -1 if elevation[0] < 0.0 else 1
    return sign * math.degrees(
        math.asin(math.fabs(elevation[0]) / math.sqrt(pow(elevation[0], 2) + pow(elevation[1], 2)))
    )


# replace method in module
astral.sun.adjust_to_obscuring_feature = adjust_to_obscuring_feature

# define location
LATITUDE = 47.8
LONGITUDE = 11.8
ELEVATION = (-10, 5000)
LOCATION = astral.Observer(LATITUDE, LONGITUDE, ELEVATION)

print(astral.sun.sun(LOCATION, date=dt.datetime.now().date())["sunrise"])
print(astral.sun.sun(LOCATION, date=dt.datetime.now().date())["sunset"])

This is what needs to work properly.

[pnbruckner]

You did not test the module as a whole calculating sunset and sunrise.

Yes, I did, I just did not share all those details here.

It comes down to how "height" (or elevation[0]) is defined. If it is defined as "height of observer relative to the obstruction", then your implementation works. However, if it is defined as "height of obstruction relative to observer", as I believe the documentation and code defines it, then it doesn't.

Both implementations "work", in that they both return a correct absolute value of the result, with the only issue being the definition of height and the needed sign that implies. I believe the implementation I suggested is a bit simpler and faster, but again, either is fine, as long as the definition of height is addressed.

[Nechoj]

You are right, it boils down to a definition. To me the line

LOCATION = astral.Observer(LATITUDE, LONGITUDE, ELEVATION)

seems to define the location of the observer, including his/her elevation. Positive elevation[0] then belongs to a higher position of the observer relative to the horizon.

Regarding the performance of the implementation, we could improve it like so:

def adjust_to_obscuring_feature(elevation: Tuple[float, float]) -> float:
    """ Calculate the number of degrees to adjust for an obscuring feature"""

    if elevation[0] == 0.0:
        return 0.0

    return math.degrees(math.asin(elevation[0]) / math.sqrt(elevation[0]*elevation[0] + elevation[1]*elevation[1]))

The handling of the sign is obsolete if asin is used, because it handles the sign already correctly. This is different when acos is used. And * is faster than pow().