skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.41k stars 212 forks source link

Confusing in indentifying the exact date day set for antitransit time in "almanac.meridian_transits()" #893

Closed khaled-hammoud closed 1 year ago

khaled-hammoud commented 1 year ago

Hi, I am using this script below to find the times for antitransit and transit but the date for transit time, exactly the returned day confuses me because I am assuming antitransit come first then transit always and not changing by changing longitude to be always in one day distance from transit time. Can you please see my script below and read comments inside when you have time to understand what I mean or is that a possible (bug) in calculation or how exactly? Thanks in advance.

from skyfield.api import load, wgs84
from skyfield import almanac
from datetime import datetime, timedelta
from pytz import timezone, utc

eph = load('de440.bsp')
earth, sun = eph['earth'], eph['sun']
ts = load.timescale()

latitude  = 52 + 20/60 + 49/3600

#A location in west poland
longitude = 15 + 43/60 + 11/3600

#A location in west netherlands offshore
#longitude = 3 + 43/60 + 11/3600   #please uncomment this longitude and comment out the previous longitude and viceversa for testing purpose

altitude = 120.0  # meters above ellipsoid WGS84
tz = timezone('Europe/Warsaw') #this is same for Poland, Germany, Belgium are the same

observatory = wgs84.latlon(latitude, longitude, altitude)

#On December 21th, 2022 at 12:00 local time
local_time = tz.localize(datetime(2022, 12, 21, 12))
utc_time = local_time.astimezone(utc) # Convert the local time to UTC
skyfield_time = ts.utc(utc_time)

observer = earth + observatory

#from
midnight = local_time.replace(hour=0)
#to
next_midnight = midnight + timedelta(days=1)

t0 = ts.from_datetime(midnight)           #on this day at 00:00:00
t1 = ts.from_datetime(next_midnight)

f = almanac.meridian_transits(eph, sun, observatory)
times, events = almanac.find_discrete(t0, t1, f)
times = times[events]

anti_transit = times[0]
transit      = times[1]

anti_transit_time = str(anti_transit.astimezone(tz))[:19]
transit_time =      str(transit.astimezone(tz))[:19]

print(anti_transit_time)
print(transit_time)

''' OUTPUT  for first place (in poland)'''
# Should it be on December 20th instead 21th ?!
# Because anti_transit occurs first then transit
# 2022-12-21 23:55:22  Anti_transit time coming AFTER transit

#it is okay
# 2022-12-21 11:55:07  transit time
#####################

''' OUTPUT for second place (offshore netherlands) '''
'''with all same input parameters but only changing LONGITUDE westward'''

# It is okay at least for me
# Because anti_transit occurs first then transit
# 2022-12-21 00:42:53   Anti_transit time coming BEFORE transit

# it is okay
# 2022-12-21 12:43:08  transit time
brandon-rhodes commented 1 year ago

Ideally the Sun's antitransit would occur at midnight, and its transit at noon, but of course time zones make our clocks deviate from that ideal. Plus, the shape of the analemma means that during some parts of the year the antitransit and transit will fall before midnight and noon, and during other parts of the year they fall after it:

https://en.wikipedia.org/wiki/Analemma

My guess is that either the time zone, or the analemma, or a combination of both, is what is moving antitransit from exactly midnight to before it (23:55:22 is, what, about 4½ minutes early)?

Try studying the whole year, maybe every day of the year in a little diagram or list of time differences, and you should see the antitransit drift one way and then the other with respect to time-zone midnight!

brandon-rhodes commented 1 year ago

(As it's now been two days, I'm going to go ahead and close this issue to keep my GitHub dashboard clear, since it looks like Skyfield is behaving correctly here; but please feel free to comment further with any other questions you might have!)

khaled-hammoud commented 1 year ago

thanks for replying, a study for one year is a good idea, so I plotted the anti_transit times for a whole year 2022, then I analyzed it and would to show you that.

Screenshot (1793) - Copy

I asked my questions on the image, exactly by the red part of curve where I suppose that it would be like this.

That is how I generated it for one year:

from skyfield.api import load, wgs84
from skyfield import almanac
from datetime import datetime, timedelta
from pytz import timezone, utc
import matplotlib.pyplot as plt

eph = load('de440.bsp')
earth, sun = eph['earth'], eph['sun']
ts = load.timescale()

latitude  = 52 + 20/60 + 49/3600
longitude = 15 + 43/60 + 11/3600
altitude = 120.0  # meters above ellipsoid WGS84
tz = timezone('Europe/Warsaw')
observatory = wgs84.latlon(latitude, longitude, altitude)

anti_transit_times = []
#transit_times      = []
day_numbers = []
for i in range(364):
    local_time = tz.localize(datetime(2022, 1, 1, 12) + timedelta(days = i))
    utc_time = local_time.astimezone(utc) # Convert the local time to UTC
    skyfield_time = ts.utc(utc_time)

    observer = earth + observatory

    #from
    midnight = local_time.replace(hour=0)
    #to
    next_midnight = midnight + timedelta(days=1)

    t0 = ts.from_datetime(midnight) #on this day at 00:00:00
    t1 = ts.from_datetime(next_midnight)

    f = almanac.meridian_transits(eph, sun, observatory)
    times, events = almanac.find_discrete(t0, t1, f)
    times = times[events]

    anti_transit = times[0]
    transit      = times[1]

    anti_transit_time = str(anti_transit.astimezone(tz))[:19]
    transit_time =      str(transit.astimezone(tz))[:19]

    ####### #Just for matplotlib purposes #################################
    day_numbers.append(datetime(2022, 1, 1) + timedelta(days = i))

    Hour,Minute = map(int,anti_transit_time[11:16].split(':'))
    anti_transit_times.append(datetime(2022, 1, 1, Hour, Minute).strftime('%H:%M'))

    #Hour,Minute = map(int,transit_time[11:16].split(':'))
    #transit_times.append(datetime(2022, 1, 1, Hour, Minute).strftime('%H:%M'))

plt.plot(day_numbers, anti_transit_times)
#plt.plot(day_numbers, transit_times)

plt.show()

#Maybe DST in Europe has to do with that?!!
#Start of DST in 2022,   2022-03-27
#End   of DST in 2022,     2022-10-30

#Here it occurs the problem betweenn these two dates
#|   anti_transit     |       transit        |
#---------------------------------------------
#|2022-10-29 00:40:52 | 2022-10-29 12:40:50  |
#|2022-10-30 23:40:45 | 2022-10-30 11:40:47  |
#---------------------------------------------
brandon-rhodes commented 1 year ago

It looks like your y-axis has a gap of one hour inside of it:

image

It's a little hard to see because the labels are so close together, but it jumps right from 00:40 to 23:40, if I'm reading it correctly. I think matplotlib must be interpreting your y axis values as simply a stack of labels, and it doesn't understand how they should be spaced because it doesn't know their real values? So it's just omitting the hour jump you want to see because of the timezone daylight savings time shift.

Similarly it looks like the plateau lines up with a label that says 01:03 or something, that's mixed in with all the times that start with 00:. So it might go away if you switch to a y-axis that matplotlib understands.

What if, instead of trying to learn how to make matplotlib understand hours-and-minutes — which might take time — you just use minutes-since-midnight (hours × 60 + minutes) as the y-axis, since this chart is just for quick debugging and to help you understand?

khaled-hammoud commented 1 year ago

I tried your method, that was good approach and plot them but it seems be confusing for me because midnight (anti_transit) come first then transit (noon).

sdsdsd
from skyfield.api import load, wgs84
from skyfield import almanac
from datetime import datetime, timedelta
from pytz import timezone, utc
import matplotlib.pyplot as plt

eph = load('de440.bsp')
earth, sun = eph['earth'], eph['sun']
ts = load.timescale()

latitude  = 52 + 20/60 + 49/3600
longitude = 15 + 43/60 + 11/3600
altitude = 120.0  # meters above ellipsoid WGS84
tz = timezone('Europe/Warsaw')
observatory = wgs84.latlon(latitude, longitude, altitude)

anti_transit_times = []
transit_times      = []
day_numbers = []
for i in range(364):
    local_time = tz.localize(datetime(2022, 1, 1, 12) + timedelta(days = i))
    utc_time = local_time.astimezone(utc) # Convert the local time to UTC
    skyfield_time = ts.utc(utc_time)

    observer = earth + observatory

    #from
    midnight = local_time.replace(hour=0)
    #to
    next_midnight = midnight + timedelta(days=1)

    t0 = ts.from_datetime(midnight) #on this day at 00:00:00
    t1 = ts.from_datetime(next_midnight)

    f = almanac.meridian_transits(eph, sun, observatory)
    times, events = almanac.find_discrete(t0, t1, f)
    times = times[events]

    anti_transit = times[0]
    transit      = times[1]

    anti_transit_time = str(anti_transit.astimezone(tz))[:19]
    transit_time =      str(transit.astimezone(tz))[:19]

    ####### #Just for matplotlib purposes #################################
    day_numbers.append(datetime(2022, 1, 1) + timedelta(days = i))

    anti_transit_time_dt_object = anti_transit.astimezone(tz)
    its_midnight = anti_transit_time_dt_object.replace(hour=0, minute=0, second=0, microsecond=0)
    anti_transit_time_min_since_midnight = (anti_transit_time_dt_object- its_midnight).total_seconds()/60.0  #minutes since midnight
    anti_transit_times.append(anti_transit_time_min_since_midnight)

    transit_time_dt_object = transit.astimezone(tz)
    its_midnight = transit_time_dt_object.replace(hour=0, minute=0, second=0, microsecond=0)
    transit_time_min_since_midnight = (transit_time_dt_object- its_midnight).total_seconds()/60.0  #minutes since midnight
    transit_times.append(transit_time_min_since_midnight)

plt.plot(day_numbers, anti_transit_times)
plt.plot(day_numbers, transit_times)

plt.show()

#Maybe DST in Europe has to do with that?!!
#Start of DST in 2022,   2022-03-27
#End   of DST in 2022,     2022-10-30

#Here it occurs the problem betweenn these two dates
#|   anti_transit     |       transit        |
#---------------------------------------------
#|2022-10-29 00:40:52 | 2022-10-29 12:40:50  |
#|2022-10-30 23:40:45 | 2022-10-30 11:40:47  |
#---------------------------------------------
brandon-rhodes commented 1 year ago

I'm glad we now have the diagram looking better! Let's now focus on the dates you asked about:

#Here it occurs the problem betweenn these two dates
#|   anti_transit     |       transit        |
#---------------------------------------------
#|2022-10-29 00:40:52 | 2022-10-29 12:40:50  |
#|2022-10-30 23:40:45 | 2022-10-30 11:40:47  |
#---------------------------------------------

Could you explain what answer you want the code to give for the date 2022-10-30? The answer it is giving looks correct to me, given that you are asking Skyfield to tell you "all the transit and anti-transit events between midnight 2022-10-30 and midnight 10-31". So I think you need to show exactly what return value you are getting for 2022-10-30, and then type out exactly what you would have expected instead for Skyfield to return.

khaled-hammoud commented 1 year ago

Let's say, I am simply a normal user of skyfield, without astronomic background etc, I just found a wonderful library like skyfield. I assume that skyfield takes into consideration DST like it is really, I would know for example when occur noon and midnight for only a given day for the location above and randomly was my date 2022-10-30 or later unti end of year 2022. I entered my input data for 2022-10-30 then skyfield return for me the noon occurs on 2022-10-30 at 11:40:47 and midnight nearly 12 hours later, then I would ask me, but midnight occurs first then noon why skyfield returns for me the noon for my given day '2022-10-30' then the midnight for the next day (end of day 2022-10-30 and start of day 2022-10-31), that would be confusing, then I start to search where is wrong with that, why starting from this date until end of year should I get as returned the midnight time of the following day of my wished and target day.

Also the last graph seems for me that it doesnt behave normal and maybe something wrong.

is there any option in skyfield to disable DST for timezone('Europe/warsaw') for example? I would redraw the graph without DST taken into consideration then I share and discuss the graph with you if both lines orange and blue remain like that or will be fart apart from each other as expected.

brandon-rhodes commented 1 year ago

So, instead of the answer that Skyfield is giving:

['2022-10-30 11:40', '2022-10-30 23:40'] [1 0]

— it sounds like you instead want this answer:

['2022-10-30 23:40', '2022-10-31 11:40'] [0 1]

If that's the case, then the problem is simply that you've forbidden Skyfield from returning that value: because you restricted Skyfield's search to the hours between one midnight and the next, it cannot return 2022-10-29 12:40, because that falls outside of the date range you are asking about.

So how about, instead of asking for events between one midnight and the next, you ask about the events between, say, 6 PM of one day, and 6 PM the next day?

    midnight = local_time.replace(hour=0) - timedelta(hours=6)

That way, whether anti-transit falls before midnight or after midnight, it will still be inside of the range of times that you are asking Skyfield to search for.

khaled-hammoud commented 1 year ago

Thank you so much, that was a good idea, I will try it later, I think you again for your time and your wonderful skyfield