我想知道是否有一个Python函数/模块，可以根据UTC时间和经度，计算午夜后的当地时间（或当地太阳时间）？ 它不需要考虑夏令时。
提前感谢。
从UTC和经度计算当地太阳时间函数
13
- Shejo284
8
5这句话的意思是：“将太阳时间与世界协调时间相差的小时数，乘以经度与π的商再乘以12，即可得到当前地点的太阳时间。” - Celada
@Celada：你有这方面的参考资料吗？我尝试了以下代码： 0 + np.arange(-180.0,190.0,15.0)/360*12
= array([-6.,-5.5,-5.,-4.5,-4.,-3.5,-3.,-2.5,-2.,-1.5,-1.,
-0.5,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.])
它只能覆盖6个小时。如果你将12小时改为24小时，你就可以得到全球的时间。如果使用pi，我不明白。你能解释一下吗？ - Shejo284
@Celada：我错过的另一件事是，经线以西的太阳时间为负数。我通过将lst中的结果进行以下处理（更像是hack）来解决这个问题：对于24小时制，让[lst < 0] += 24。我相信一个更懂数学的人可以用一行代码解决这个问题。 - Shejo284
1@J.F.Sebastian，我不知道这个问题如此复杂。在Wikipedia上进行了一些挖掘，显示当地时间在一年中至少会漂移16分钟的示例，尽管它在长期内平均化。 - Mark Ransom
@Celada，这并不是显而易见的 - 这就是我为什么说的原因。我曾经看到人们混淆度和弧度，而且他们并不知道这一点。无论电脑使用什么内部表示方式，可读的经度测量值总是以度为单位给出。附言：J.F. Sebastian的答案将你的简单解决方案与更准确的解决方案进行了比较，他使用的库似乎使用弧度。 - Mark Ransom
显示剩余3条评论
5个回答
14
使用ephem的sidereal_time()方法：
import ephem # pip install pyephem (on Python 2)
# pip install ephem (on Python 3)
def solartime(observer, sun=ephem.Sun()):
sun.compute(observer)
# sidereal time == ra (right ascension) is the highest point (noon)
hour_angle = observer.sidereal_time() - sun.ra
return ephem.hours(hour_angle + ephem.hours('12:00')).norm # norm for 24h
注意：
ephem.hours是一个浮点数，表示弧度角，并将其转换为/从字符串“hh：mm：ss.ff”。为了比较，这里是“UTC +经度”公式：
import math
from datetime import timedelta
def ul_time(observer):
utc_dt = observer.date.datetime()
longitude = observer.long
return utc_dt + timedelta(hours=longitude / math.pi * 12)
例子
from datetime import datetime
# "solar time" for some other cities
for name in ['Los Angeles', 'New York', 'London',
'Paris', 'Moscow', 'Beijing', 'Tokyo']:
city = ephem.city(name)
print("%-11s %11s %s" % (name, solartime(city),
ul_time(city).strftime('%T')))
# set date, longitude manually
o = ephem.Observer()
o.date = datetime(2012, 4, 15, 1, 0, 2) # some utc time
o.long = '00:00:00.0' # longitude (you could also use a float (radians) here)
print("%s %s" % (solartime(o), ul_time(o).strftime('%T')))
输出
Los Angeles 14:59:34.11 14:44:30
New York 17:56:31.27 17:41:27
London 22:52:02.04 22:36:58
Paris 23:01:56.56 22:46:53
Moscow 1:23:00.33 01:07:57
Beijing 6:38:09.53 06:23:06
Tokyo 8:11:17.81 07:56:15
1:00:00.10 01:00:01
- jfs
2
2这个答案很老了，所以
ephem 可能已经改变了，但是现在 solartime(o) 返回一个 float。 - duff182@duff18：它返回
ephem.Angle。它是一个 float 实例，如答案中所述，以 "hh:mm:ss.f" 的形式打印。尝试使用 str(solartime(o)) 查看它不仅仅是一个浮点数。 - jfs8
如果你想更简单一些，你可以使用NOAA的低精度公式：
#!/usr/local/bin/python
import sys
from datetime import datetime, time, timedelta
from math import pi, cos, sin
def solar_time(dt, longit):
return ha
def main():
if len(sys.argv) != 4:
print 'Usage: hour_angle.py [YYYY/MM/DD] [HH:MM:SS] [longitude]'
sys.exit()
else:
dt = datetime.strptime(sys.argv[1] + ' ' + sys.argv[2], '%Y/%m/%d %H:%M:%S')
longit = float(sys.argv[3])
gamma = 2 * pi / 365 * (dt.timetuple().tm_yday - 1 + float(dt.hour - 12) / 24)
eqtime = 229.18 * (0.000075 + 0.001868 * cos(gamma) - 0.032077 * sin(gamma) \
- 0.014615 * cos(2 * gamma) - 0.040849 * sin(2 * gamma))
decl = 0.006918 - 0.399912 * cos(gamma) + 0.070257 * sin(gamma) \
- 0.006758 * cos(2 * gamma) + 0.000907 * sin(2 * gamma) \
- 0.002697 * cos(3 * gamma) + 0.00148 * sin(3 * gamma)
time_offset = eqtime + 4 * longit
tst = dt.hour * 60 + dt.minute + dt.second / 60 + time_offset
solar_time = datetime.combine(dt.date(), time(0)) + timedelta(minutes=tst)
print solar_time
if __name__ == '__main__':
main()
- mattexx
6
1奇怪，不知道那是什么意思。是的，这对于连续的美国工作得很好。但对于阿拉斯加来说就不太行了。基本上适用于大多数纬度较小的东西。 - mattexx
1已转换为Python3 https://gist.github.com/mnpenner/8b6e7cceb930402d3f37ba9c59fa1d26。谢谢！ - mpen
2decl似乎已被定义但从未使用？ - Sebastian
1链接的PDF文件中提到，在计算闰年的伽马值时应使用366，因此您需要...
if dt.year % 4 == 0 and (dt.year % 100 != 0 or dt.year % 400 == 0):
yrlen = 366
else:
yrlen = 365
gamma = 2 * pi / yrlen * (dt.timetuple().tm_yday - 1 + float(dt.hour - 12) / 24) - Darius4
我查阅了Jean Meeus的《天文算法》。我认为你可能需要询问本地时角，可以用时间（0-24小时）、度数（0-360）或弧度（0-2pi）表示。
我猜你可以使用ephem库来实现这一功能。但是为了好玩，这里还有一些Python代码：
#!/usr/local/bin/python
import sys
from datetime import datetime, time, timedelta
from math import pi, sin, cos, atan2, asin
# helpful constant
DEG_TO_RAD = pi / 180
# hardcode difference between Dynamical Time and Universal Time
# delta_T = TD - UT
# This comes from IERS Bulletin A
# ftp://maia.usno.navy.mil/ser7/ser7.dat
DELTA = 35.0
def coords(yr, mon, day):
# @input year (int)
# @input month (int)
# @input day (float)
# @output right ascention, in radians (float)
# @output declination, in radians (float)
# get julian day (AA ch7)
day += DELTA / 60 / 60 / 24 # use dynamical time
if mon <= 2:
yr -= 1
mon += 12
a = yr / 100
b = 2 - a + a / 4
jd = int(365.25 * (yr + 4716)) + int(30.6 * (mon + 1)) + day + b - 1524.5
# get sidereal time at greenwich (AA ch12)
t = (jd - 2451545.0) / 36525
# Calculate mean equinox of date (degrees)
l = 280.46646 + 36000.76983 * t + 0.0003032 * t**2
while (l > 360):
l -= 360
while (l < 0):
l += 360
# Calculate mean anomoly of sun (degrees)
m = 357.52911 + 35999.05029 * t - 0.0001537 * t**2
# Calculate eccentricity of Earth's orbit
e = 0.016708634 - 0.000042037 * t - 0.0000001267 * t**2
# Calculate sun's equation of center (degrees)
c = (1.914602 - 0.004817 * t - .000014 * t**2) * sin(m * DEG_TO_RAD) \
+ (0.019993 - .000101 * t) * sin(2 * m * DEG_TO_RAD) \
+ 0.000289 * sin(3 * m * DEG_TO_RAD)
# Calculate the sun's radius vector (AU)
o = l + c # sun's true longitude (degrees)
v = m + c # sun's true anomoly (degrees)
r = (1.000001018 * (1 - e**2)) / (1 + e * cos(v * DEG_TO_RAD))
# Calculate right ascension & declination
seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * 0.001813))
e0 = 23 + (26 + (seconds / 60)) / 60
ra = atan2(cos(e0 * DEG_TO_RAD) * sin(o * DEG_TO_RAD), cos(o * DEG_TO_RAD)) # (radians)
decl = asin(sin(e0 * DEG_TO_RAD) * sin(o * DEG_TO_RAD)) # (radians)
return ra, decl
def hour_angle(dt, longit):
# @input UTC time (datetime)
# @input longitude (float, negative west of Greenwich)
# @output hour angle, in degrees (float)
# get gregorian time including fractional day
y = dt.year
m = dt.month
d = dt.day + ((dt.second / 60.0 + dt.minute) / 60 + dt.hour) / 24.0
# get right ascention
ra, _ = coords(y, m, d)
# get julian day (AA ch7)
if m <= 2:
y -= 1
m += 12
a = y / 100
b = 2 - a + a / 4
jd = int(365.25 * (y + 4716)) + int(30.6 * (m + 1)) + d + b - 1524.5
# get sidereal time at greenwich (AA ch12)
t = (jd - 2451545.0) / 36525
theta = 280.46061837 + 360.98564736629 * (jd - 2451545) \
+ .000387933 * t**2 - t**3 / 38710000
# hour angle (AA ch13)
ha = (theta + longit - ra / DEG_TO_RAD) % 360
return ha
def main():
if len(sys.argv) != 4:
print 'Usage: hour_angle.py [YYYY/MM/DD] [HH:MM:SS] [longitude]'
sys.exit()
else:
dt = datetime.strptime(sys.argv[1] + ' ' + sys.argv[2], '%Y/%m/%d %H:%M:%S')
longit = float(sys.argv[3])
ha = hour_angle(dt, longit)
# convert hour angle to timedelta from noon
days = ha / 360
if days > 0.5:
days -= 0.5
td = timedelta(days=days)
# make solar time
solar_time = datetime.combine(dt.date(), time(12)) + td
print solar_time
if __name__ == '__main__':
main()
- mattexx
3
谢谢您的回复，但您的方向有些偏离。我想要的是一个简单的公式，可以根据经纬度和UTC时间计算出当地的太阳时间。 - Shejo284
那么我的答案（以度为单位的时间角）不是你要找的，还是解决方案太复杂了？ - mattexx
1我更新了代码，现在可以生成太阳时间作为日期时间，以防这是你问题的一部分。我会尝试在另一个答案中使用ephem来重现这个问题，因为我猜想这可能只是比你要寻找的更多的代码行。 - mattexx
3
对于一个非常简单的功能和非常近似的本地时间：
时间变化范围从-12小时到+12小时，经度范围从-180到180。然后：
import datetime as dt
def localTimeApprox(myDateTime, longitude):
"""Returns local hour approximation"""
return myDateTime+dt.timedelta(hours=(longitude*12/180))
示例调用：localTimeApprox(dt.datetime(2014, 7, 9, 20, 00, 00), -75)
返回：datetime.datetime(2014, 7, 9, 15, 0)
- Le Droid
4
Droid：优雅！谢谢！ - Shejo284
@Shejo284：-1. 这很简单（与我答案中的
ul_time()相同，只是它接受日期时间、经度（以度为单位）而不是提供相同信息的观察者），但是这是错误的。它不显示您的当地太阳时间（就像简单的日晷一样）。将其结果与solartime()函数进行比较。我的答案末尾的输出显示全年约有15分钟的差异。 - jfs1@J.F.Sebastian 我应该编辑我的答案，在“非常非常近似”的“非常”后再加一个吗？只需要15分钟？这比我想象的要好，谢谢。如果你觉得安装额外的库并编写超过10行代码比一行代码更简单，那是你的意见，但对我来说不是一个选项，因为我不能没有一个真正好的理由在每个服务器上都安装它。至于给聪明的解决方案罚款，嗯...有趣。 - Le Droid
你是指我回答中复制
ul_time() 的“聪明”解决方案吗？ - jfs2
再次尝试，使用ephem。我将纬度和海拔作为参数包含在内，但当然不需要。你可以为你的目的调用它们为0。
#!/usr/local/bin/python
import sys
from datetime import datetime, time, timedelta
import ephem
def hour_angle(dt, longit, latit, elev):
obs = ephem.Observer()
obs.date = dt.strftime('%Y/%m/%d %H:%M:%S')
obs.lon = longit
obs.lat = latit
obs.elevation = elev
sun = ephem.Sun()
sun.compute(obs)
# get right ascention
ra = ephem.degrees(sun.g_ra) - 2 * ephem.pi
# get sidereal time at greenwich (AA ch12)
jd = ephem.julian_date(dt)
t = (jd - 2451545.0) / 36525
theta = 280.46061837 + 360.98564736629 * (jd - 2451545) \
+ .000387933 * t**2 - t**3 / 38710000
# hour angle (AA ch13)
ha = (theta + longit - ra * 180 / ephem.pi) % 360
return ha
def main():
if len(sys.argv) != 6:
print 'Usage: hour_angle.py [YYYY/MM/DD] [HH:MM:SS] [longitude] [latitude] [elev]'
sys.exit()
else:
dt = datetime.strptime(sys.argv[1] + ' ' + sys.argv[2], '%Y/%m/%d %H:%M:%S')
longit = float(sys.argv[3])
latit = float(sys.argv[4])
elev = float(sys.argv[5])
# get hour angle
ha = hour_angle(dt, longit, latit, elev)
# convert hour angle to timedelta from noon
days = ha / 360
if days > 0.5:
days -= 0.5
td = timedelta(days=days)
# make solar time
solar_time = datetime.combine(dt.date(), time(12)) + td
print solar_time
if __name__ == '__main__':
main()
- mattexx
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接