search

skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.38k stars 208 forks source link

Document that Python datetime cannot represent Julian calendar leap days like February 29 of the year 1000 #957

Open godbolerr opened 2 months ago

godbolerr commented 2 months ago

I am calculating sunrise for my location from year 1000 onwards.

Skyfield is unable to compute the same and I got the following error. It works for dates before and after this value 1000-3-1.

Here is the stacktrace of the error

Traceback (most recent call last): File "E:\Astronomy\AstroExhibits\SkyfieldTutorial\LeapSecondInSkyfield.py", line 36, in sunriseLocal = st[0].astimezone(istZone) ^^^^^^^^^^^^^^^^^^^^^^^^^ File "F:\Python\install\Lib\site-packages\skyfield\timelib.py", line 454, in astimezone dt, leap_second = self.astimezone_and_leap_second(tz) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "F:\Python\install\Lib\site-packages\skyfield\timelib.py", line 480, in astimezone_and_leap_second dt, leap_second = self.utc_datetime_and_leap_second() ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "F:\Python\install\Lib\site-packages\skyfield\timelib.py", line 544, in utc_datetime_and_leap_second dt = datetime(year, month, day, hour, minute, second, micro, utc) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ValueError: day is out of range for month

Here is the sample program to reproduce the error.

from datetime import datetime from datetime import timedelta import datetime

from pytz import timezone import pytz
from skyfield import api from skyfield import almanac from skyfield.api import GREGORIAN_START from skyfield.api import N, E, wgs84, load

ts = api.load.timescale() istZone = timezone('Asia/Kolkata') ts.julian_calendar_cutoff = GREGORIAN_START eph = load('de431t.bsp') lat=18.5204 long=73.8567 sun, moon, earth = eph['sun'], eph['moon'], eph['earth'] pune = wgs84.latlon(lat N, long E) puneObserver = eph['Earth'] + pune istZone = timezone('Asia/Kolkata')

refDate = datetime.datetime(1000,3,1,0,0,0) midnight = refDate.replace(tzinfo=istZone) next_midnight = midnight + datetime.timedelta(days=1)

t0 = ts.from_datetime(midnight) t1 = ts.from_datetime(next_midnight)

st, _ = almanac.find_risings(puneObserver, sun, t0, t1)

sunriseLocal = st[0].astimezone(istZone)

print(sunriseLocal)

brandon-rhodes commented 2 months ago

You need to move away from the Python datetime to see that date. Let's try the robust JPL way of printing the date:

print(st[0].utc_jpl())

The date is:

A.D. 1000-Feb-29 01:48:10.4650 UTC

Python only understands the Gregorian calendar, in which there is no leap day in the year 1000 per the rules. But you are using the historic Julian calendar for old dates thanks to ts.julian_calendar_cutoff = GREGORIAN_START, and in the Julian calendar — whose leap year rules are simpler — there was indeed a leap day in the year 1000.

So you'll have to avoid Python's datetime.

Let's keep this issue open until I have time to write up docs about this. Maybe I should even add detection for this error, so that folks get a message that tells them exactly what went wrong.

Thanks for the report!