search

pingswept / pysolar

Pysolar is a collection of Python libraries for simulating the irradiation of any point on earth by the sun. It includes code for extremely precise ephemeris calculations.
http://pysolar.org
GNU General Public License v3.0
372 stars 125 forks source link

util.py/get_sunrise_sunset_transit() uses the wrong day? #144

Open rm-minus-r-star opened 3 years ago

rm-minus-r-star commented 3 years ago

I ran into a problem using get_sunrise_sunset_transit() in util.py. "same_day" goes to the beginning of the day for the input time, not for the day that's relevant at the specified longitude. This might be desired behaviour with respect to other functions within the package that depend on this. But for my purposes, I work in UTC to avoid ambiguity over DST, but determine sunrise/sunset times and solar angles for a timezone that's +10 from UTC. If I use a UTC time of, say, 2021-12-31 23:00:00+00:00, Sunrise should be more like 2021-12-31 19:03:22+00:00, not 2021-12-30 19:03:22+00:00 2021, 12, 30, 19, 3, 22+00:00

For now I've replaced lines 159-162 with the following:

    time_diff = int(round(longitude_deg / 15., 0) - utc_offset / 3600)
    local_day = when + timedelta(hours = time_diff)
    local_day = datetime(year = local_day.year, month = local_day.month, day = local_day.day, tzinfo = local_day.tzinfo)
    print("local_day", local_day)

    transit_time = local_day + timedelta(hours = time_diff + TON)
    transit_time = (transit_time - timedelta(hours = time_diff)).replace(tzinfo=when.tzinfo)
    sunrise_time = transit_time - timedelta(hours = ha)
    sunset_time = transit_time + timedelta(hours = ha)

I think this is right, but I'm not sure of the physics ... I'm not sure whether the earlier day = when.timetuple().tm_yday # Day of the year needs to be similarly adjusted according to the longitude or not.

pingswept commented 3 years ago

Yeah, it seems like you've identified a problem. get_sunrise_sunset_transit() should check the local time to determine the day, and then determine the sunrise, sunset, and transit times from that. I would welcome a patch that did that.

Luigolas commented 2 years ago

We have recently updated out pysolar version from 0.9 to 0.10 and suddenly few tests are broken. I'm struggling to understand if this is the correct behaviour but when doing this:

import pytz
from datetime import datetime
import pysolar.util

day = datetime(2018, 6, 21, 0, 0, tzinfo=pytz.timezone('Europe/Zurich'))
lat, lon = 47.5642388562567, 7.636568666258004
pysolar.util.get_sunrise_sunset(latitude_deg=lat, longitude_deg=lon, when=day)

>> (datetime.datetime(2018, 6, 20, 4, 39, 28, 511700, tzinfo=<DstTzInfo 'Europe/Zurich' LMT+0:34:00 STD>), datetime.datetime(2018, 6, 20, 20, 38, 8, 861432, tzinfo=<DstTzInfo 'Europe/Zurich' LMT+0:34:00 STD>))

As you can see I get the previous day, not the 21. Tested also setting the hour at 12, but then results are shifted:

pysolar.util.get_sunrise_sunset(latitude_deg=lat, longitude_deg=lon, when=day.replace(hour=12))
>> (datetime.datetime(2018, 6, 21, 4, 39, 38, 510017, tzinfo=<DstTzInfo 'Europe/Zurich' LMT+0:34:00 STD>), datetime.datetime(2018, 6, 21, 20, 38, 24, 952483, tzinfo=<DstTzInfo 'Europe/Zurich' LMT+0:34:00 STD>))

In order to make it work (as in 0.9) I have to:

  1. Make sure the datetime object is in the middle of the day.
  2. Convert to UTC before passing it to the _get_sunrisesunset method
  3. Convert results to local timezone

Am I doing something wrong here?

bhemmer commented 2 years ago

I have the same issue, any news on that?

jhaskinsPhD commented 1 year ago

Fixed this error if you pass "when" in UTC (somehow it's still breaking for utcoffset!=0, and I don't have time to fix that too). It was calculating sunrise/sunset using transit_time in UTC, not in local time. Because ha > transit_time.hour, often if transit_time is in UTC, and you didn't call the function "in the middle of the day", that was causing the sunrise time to end up on the previous day. This can be fixed by making sure you calculate sunrise/sunset time from transit_time in local time, then converting those back to the time zone they were passed in. This gives me back good sunrise/sunset times in UTC that make sense if I pass "when" in UTC.

Lines 167-181

# Get the local day. 
local_day = datetime(year = local_time.year, month = local_time.month, day = local_time.day, tzinfo = local_time.tzinfo)

# Define transit time in Local time (tzinfo says UTC even tho its in local time
# because tzinfo for local_time is never defined in this function, so its wrong, but doesn't matter here. 
# We know transit time is in local time because it is, and should be around 12. 
transit_time = local_day + timedelta(hours = time_diff + TON) 

# Now get sunrise & sunset in local time! 
sunrise_time_local = transit_time - timedelta(hours = ha)    
sunset_time_local = transit_time + timedelta(hours = ha)

# Finally, convert it back to the native time zone it was passed in! 
# Note: time_diff is always positive, while utc_offset can actually be negative- that's why the signs are diff! 
sunrise_time=sunrise_time_local +timedelta(hours=time_diff)-timedelta(seconds=utc_offset)
sunset_time=sunset_time_local +timedelta(hours=time_diff)-timedelta(seconds=utc_offset)
wolfgang-302 commented 1 year ago

i dont know, iff i found the same error. But

from datetime import datetime
from zoneinfo import ZoneInfo
zeit = datetime(2023,5,1,tzinfo=ZoneInfo('Europe/Berlin'))
util.get_transit_time(51,7,zeit).time()

leads to datetime.time(15, 29, 10, 541377) which has a shift of 2 hours.

zeit = datetime(2023,5,1,tzinfo=ZoneInfo('utc'))
util.get_transit_time(51,7,zeit).astimezone(ZoneInfo('Europe/Berlin')).time()

leads to the correct value datetime.time(13, 29, 2, 938833)

I had a look at the sources but i dont really understand what is calculated there. So i came up with my own function get_transit(longitude,d) which is defined as

from datetime import datetime, timedelta, time
from zoneinfo import ZoneInfo
from pysolar import solar

def get_transit(longitude,d):
    ''' return the time of transit of the sun 
        for a given longitude and a given datetime d
        If d is not timezoneaware, the local
        timezone is used
    '''

    tz = d.tzinfo if d.tzinfo is not None else d.astimezone().tzinfo
    eot = solar.equation_of_time(d.utctimetuple().tm_yday)
    ret_val = d.replace(hour=12,tzinfo=ZoneInfo('utc'))
    ret_val -= timedelta(minutes=eot) # correction wrt equation of time
    ret_val -= timedelta(hours=longitude/15) # correction wrt longitude

    return ret_val.astimezone(tz)

I did not test this code for every timezone and for north / south hemisphere. But maybe it helps to analyse the problem with util.get_transit_time().

The result has the same degree of accuracy as util.get_transit_time().

zeit = datetime(2023,5,1,tzinfo=ZoneInfo('Europe/Berlin'))
tt = get_transit(7,zeit).time()
tt

which leads to datetime.time(13, 29, 11, 993236)

SelmaUrban commented 1 year ago

Hi, you can change

    sunrise_time = transit_time - timedelta(hours = ha)
    sunset_time = transit_time + timedelta(hours = ha)

to

    sunrise_time = (transit_time - timedelta(hours = ha)) - timedelta(seconds = utc_offset)
    sunset_time = (transit_time + timedelta(hours = ha)) - timedelta(seconds = utc_offset)

That fixed the problem for me.

Kuenlun commented 2 days ago

The function get_sunrise_sunset_transit worked well previous to the fix added by bf02e56d7412ce9910ba084ef538f8276acc36dc. To demonstrate it I've created 3 python test cases and an additional C test with the original algorithm. The python tests show:

  1. The output of the function now (28c77b0e91c8674e3d755f382cf5d112763d8719 commited in Aug 05 2024, previous to the pull request I just opened).
  2. The output of the function after the fix I made almost 4 years ago (3de002288dd8971ec11587c05694776709e8b775 commited in Nov 29 2020)
  3. The output of the function previous to my first fix (2d4b89cd6243bca051680f3ff50b8b311bf51ccf commited in Oct 18 2020)

The correct output, at least compared to the original C algorithm is the number 2.

For more information about the fix I made 4 years ago and why did I propose it see #137.

Python test

from datetime import datetime, time, date
from pysolar.util import get_transit_time

# Get the UTC offset
utc_offset = datetime.now().astimezone().utcoffset()
print(f'UTC offset is {utc_offset}')

# Date to use as an example
day = date(year=2003, month=10, day=17)

# Iterate over 24 hours and calculate noon at each moment
for hour in range(24):
    when = datetime.combine(day, time(hour=hour, minute=30, second=30)).astimezone()
    transit = get_transit_time(39.742476, -105.1786, when)
    print(f'Transit calculated at {when} is {transit}')

Output at last commit (28c77b0e91c8674e3d755f382cf5d112763d8719) commited in Aug 05 2024:

UTC offset is 2:00:00
Transit calculated at 2003-10-17 00:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 01:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 02:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 03:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 04:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 05:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 06:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 07:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 08:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 09:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 10:30:30+02:00 is 2003-10-16 22:46:14.151396+02:00
Transit calculated at 2003-10-17 11:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 12:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 13:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 14:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 15:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 16:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 17:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 18:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 19:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 20:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 21:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 22:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00
Transit calculated at 2003-10-17 23:30:30+02:00 is 2003-10-17 22:46:01.552022+02:00

Output just after my fix (3de002288dd8971ec11587c05694776709e8b775) commited in Nov 29 2020:

UTC offset is 2:00:00
Transit calculated at 2003-10-17 00:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 01:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 02:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 03:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 04:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 05:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 06:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 07:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 08:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 09:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 10:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 11:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 12:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 13:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 14:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 15:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 16:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 17:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 18:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 19:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 20:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 21:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 22:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 23:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00

Output previous to my first fix (2d4b89cd6243bca051680f3ff50b8b311bf51ccf) commited in Oct 18 2020:

UTC offset is 2:00:00
Transit calculated at 2003-10-17 00:30:30+02:00 is 2003-10-17 20:46:14.151396+02:00
Transit calculated at 2003-10-17 01:30:30+02:00 is 2003-10-17 20:46:14.151396+02:00
Transit calculated at 2003-10-17 02:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 03:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 04:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 05:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 06:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 07:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 08:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 09:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 10:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 11:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 12:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 13:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 14:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 15:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 16:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 17:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 18:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 19:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 20:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 21:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 22:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00
Transit calculated at 2003-10-17 23:30:30+02:00 is 2003-10-17 20:46:01.552022+02:00

C test

int main (int argc, char *argv[])
{
    spa_data spa;  //declare the SPA structure
    int result;
    float min, sec;

    //enter required input values into SPA structure
    for (int i{ 0 }; i < 24; ++i)
    {
        spa.year = 2003;
        spa.month = 10;
        spa.day = 17;
        spa.hour = i;
        spa.minute = 30;
        spa.second = 30;
        spa.timezone = +2.0;
        spa.delta_ut1 = 0;
        spa.delta_t = 67;
        spa.longitude = -105.1786;
        spa.latitude = 39.742476;
        spa.elevation = 1830.14;
        spa.pressure = 820;
        spa.temperature = 11;
        spa.slope = 30;
        spa.azm_rotation = -10;
        spa.atmos_refract = 0.5667;
        spa.function = SPA_ALL;

        //call the SPA calculate function and pass the SPA structure
        result = spa_calculate(&spa);
        if (result == 0)  //check for SPA errors
        {
            //display the results inside the SPA structure
            min = 60.0 * (spa.suntransit - (int)(spa.suntransit));
            sec = 60.0 * (min - (int)min);
            printf("Transit:       %02d:%02d:%02d Local Time\n", (int)(spa.suntransit), (int)min, (int)sec);
        } else printf("SPA Error Code: %d\n", result);
    }

    return 0;
}

Output from C test:

Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time
Transit:       20:46:04 Local Time