日出/日落计算

7

我正在尝试使用Python根据下面提供的链接计算日落/日出时间。

我的Excel和Python计算结果与实际值不匹配。您有什么想法,我可能做错了什么吗?

我的Excel表格可以在以下位置找到... http://transpotools.com/sun_time.xls

# Created on 2010-03-28

# @author: dassouki
# @source: [http://williams.best.vwh.net/sunrise_sunset_algorithm.htm][2]
# @summary: this is based on the Nautical Almanac Office, United States Naval
# Observatory.

import math, sys

class TimeOfDay(object):

  def calculate_time(self, in_day, in_month, in_year,
                     lat, long, is_rise, utc_time_zone):

    # is_rise is a bool when it's true it indicates rise,
    # and if it's false it indicates setting time

    #set Zenith
    zenith = 96
    # offical      = 90 degrees 50'
    # civil        = 96 degrees
    # nautical     = 102 degrees
    # astronomical = 108 degrees


    #1- calculate the day of year
    n1 = math.floor( 275 * in_month / 9 )
    n2 = math.floor( ( in_month + 9 ) / 12 )
    n3 = ( 1 + math.floor( in_year - 4 * math.floor( in_year / 4 ) + 2 ) / 3 )

    new_day = n1 - ( n2 * n3 ) + in_day - 30

    print "new_day ", new_day

    #2- calculate rising / setting time
    if is_rise:
      rise_or_set_time = new_day + ( ( 6 - ( long / 15 ) ) / 24 )
    else:
      rise_or_set_time = new_day + ( ( 18 - ( long/ 15 ) ) / 24 )

    print "rise / set", rise_or_set_time

    #3- calculate sun mean anamoly
    sun_mean_anomaly = ( 0.9856 * rise_or_set_time ) - 3.289
    print "sun mean anomaly", sun_mean_anomaly

    #4 calculate true longitude
    true_long = ( sun_mean_anomaly +
                  ( 1.916 * math.sin( math.radians( sun_mean_anomaly ) ) ) +
                  ( 0.020 * math.sin(  2 * math.radians( sun_mean_anomaly ) ) ) +
                  282.634 )
    print "true long ", true_long

    # make sure true_long is within 0, 360
    if true_long < 0:
      true_long = true_long + 360
    elif true_long > 360:
      true_long = true_long - 360
    else:
      true_long

    print "true long (360 if) ", true_long

    #5 calculate s_r_a (sun_right_ascenstion)
    s_r_a = math.degrees( math.atan( 0.91764 * math.tan( math.radians( true_long ) ) ) )
    print "s_r_a is ", s_r_a

    #make sure it's between 0 and 360
    if s_r_a < 0:
      s_r_a = s_r_a + 360
    elif true_long > 360:
      s_r_a = s_r_a - 360
    else:
      s_r_a

    print "s_r_a (modified) is ", s_r_a

    # s_r_a has to be in the same Quadrant as true_long
    true_long_quad = ( math.floor( true_long / 90 ) ) * 90
    s_r_a_quad = ( math.floor( s_r_a / 90 ) ) * 90
    s_r_a = s_r_a + ( true_long_quad - s_r_a_quad )

    print "s_r_a (quadrant) is ", s_r_a

    # convert s_r_a to hours
    s_r_a = s_r_a / 15

    print "s_r_a (to hours) is ", s_r_a


    #6- calculate sun diclanation in terms of cos and sin
    sin_declanation = 0.39782 * math.sin( math.radians ( true_long ) )
    cos_declanation = math.cos( math.asin( sin_declanation ) )

    print " sin/cos declanations ", sin_declanation, ", ", cos_declanation

    # sun local hour
    cos_hour = ( math.cos( math.radians( zenith ) ) -
                 ( sin_declanation * math.sin( math.radians ( lat ) ) ) /
                 ( cos_declanation * math.cos( math.radians ( lat ) ) ) )

    print "cos_hour ", cos_hour

    # extreme north / south
    if cos_hour > 1:
      print "Sun Never Rises at this location on this date, exiting"
      # sys.exit()
    elif cos_hour < -1:
      print "Sun Never Sets at this location on this date, exiting"
      # sys.exit()

    print "cos_hour (2)", cos_hour


    #7- sun/set local time calculations
    if is_rise:
      sun_local_hour =  ( 360 - math.degrees(math.acos( cos_hour ) ) ) / 15
    else:
      sun_local_hour = math.degrees( math.acos( cos_hour ) ) / 15


    print "sun local hour ", sun_local_hour

    sun_event_time = sun_local_hour + s_r_a - ( 0.06571 *
                                                rise_or_set_time ) - 6.622

    print "sun event time ", sun_event_time

    #final result
    time_in_utc =  sun_event_time - ( long / 15 ) + utc_time_zone

    return time_in_utc



#test through main
def main():
  print "Time of day App "
  # test: fredericton, NB
  # answer: 7:34 am
  long = 66.6
  lat = -45.9
  utc_time = -4
  d = 3
  m = 3
  y = 2010
  is_rise = True

  tod = TimeOfDay()
  print "TOD is ", tod.calculate_time(d, m, y, lat, long, is_rise, utc_time)


if __name__ == "__main__":
  main()

不匹配实际值。在哪里?您打印了许多中间结果。哪些不匹配? - S.Lott
time_in_utc这个名称有误导性。显然,该函数的意图是返回本地(针对给定时区)时间。顺便说一下,电子表格显示日出时间为7:02(而不是您代码注释中的7:34 am)。 - jfs
3个回答

10

您可以使用ephem Python 模块:

#!/usr/bin/env python
import datetime
import ephem # to install, type$ pip install pyephem

def calculate_time(d, m, y, lat, long, is_rise, utc_time):
    o = ephem.Observer()
    o.lat, o.long, o.date = lat, long, datetime.date(y, m, d)
    sun = ephem.Sun(o)
    next_event = o.next_rising if is_rise else o.next_setting
    return ephem.Date(next_event(sun, start=o.date) + utc_time*ephem.hour
                      ).datetime().strftime('%H:%M')

例子:

for town, kwarg in { "Fredericton": dict(d=3, m=3, y=2010,
                                         lat='45.959045', long='-66.640509',
                                         is_rise=True,
                                         utc_time=20),

                     "Beijing": dict(d=29, m=3, y=2010,
                                     lat='39:55', long='116:23',
                                     is_rise=True,
                                     utc_time=+8),

                     "Berlin": dict(d=4, m=4, y=2010,
                                    lat='52:30:2', long='13:23:56',
                                    is_rise=False,
                                    utc_time=+2) ,

                     "Moscow": dict(d=4, m=4, y=2010,
                                    lat='55.753975', long='37.625427',
                                    is_rise=True,
                                    utc_time=4) }.items():
    print town, calculate_time(**kwarg)

输出:

Beijing 06:02
Berlin 19:45
Moscow 06:53
Fredericton 07:01

非常感谢,我将使用您的代码来实现我的项目。有时候Python就像苹果一样,无所不能。Python中总会有一个库可以解决问题。 - dassouki

3
为什么要调用“radians”和“degrees”?我以为输入数据已经是十进制度数了。
如果我执行以下操作,将所有对弧度和角度的调用删除,将纬度/经度更正为45.9和-66.6,并将time_in_utc更正为0到24之间的时间,则会得到7:37 am的结果。

编辑:
正如J. F. Sebastian所指出的那样,根据问题中链接的电子表格和使用ephem的Observer类提供的答案,在这个位置的日出时间在07:01-07:02左右。

一旦我得到了一个大致正确的数字(实现中评论中的07:34),我就停止了对dassouki的美国海军天文台算法实现中错误的寻找。

研究发现,该算法进行了一些简化,并且有关“日出”构成的变化是有差异的,其中一些在此处讨论。然而,根据我最近在这个问题上学到的知识,这些变化只应导致日出时间相差几分钟,而不是半个多小时。


根据电子表格,日出时间应该是7:02,而根据Python模块ephem,应该是7:01(而不是您之前回答的7:37)。请参考https://dev59.com/y0zSa4cB1Zd3GeqPqNyK#2539597。 - jfs
糟糕。在接近问题中指定的答案后,我停止了寻找。 - MattH
这里的惯例是什么?我应该删除我的答案吗?从某种角度来看,我已经成功回答了“我可能做错了什么?”的问题,但显然从另一个角度来看,计算出的结果并不正确。 - MattH
@MattH:评论是可见的,所以现在这样就可以了。你可以编辑你的答案指出这个差异。 - jfs

1
我怀疑这与实际上没有执行浮点除法有关。在Python中,如果a和b都是整数,则a / b也是一个整数:
   $ python
   >>> 1 / 2
   0

你的选择可以是将你的参数强制转换为浮点数(即,将a/b改为float(a) / b),或者确保“/”以Python 3K的方式正常工作:

   $ python
   >>> from __future__ import division
   >>> 1 / 2
   0.5

所以,如果你把这个导入语句放在文件的顶部,它可能会解决你的问题。现在 / 操作符将始终产生一个浮点数,如果你想要旧的行为,可以使用 // 操作符。


我的值都不是整数,但是你的建议并没有解决问题。 - dassouki
嗯,我就怕这个;-) 跟踪任何错误的唯一方法是真正地逐步执行代码,不幸的是。你确定有正确的答案可以进行比较吗? - rlotun

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接