astronomy - 使用 PyEphem 在世界地图上获取太阳的纬度和经度

Question

我正在尝试确定太阳、月亮和火星的纬度和经度。我需要相对于地球赤道和本初子午线的结果，以便产生类似于 this map 的结果。 .

我相信this question的作者也是如此。想要的，但是那里的答案对我来说并不重要(与第一个链接中的值相比)。

预期结果，从page linked to earlier获得:

On Thursday, 1 January 2015, 00:00:00 UTC the Sun is at its zenith at Latitude: 23° 02' South, Longitude: 179° 29' West

>>> import ephem; from math import degrees
>>> b = ephem.Sun(epoch='date'); b.compute('2015/1/1 00:00:00')
>>> print("{},{}".format(degrees(b.dec), degrees(b.ra)))
-23.040580418272267,281.12827017399906

因此，纬度/赤纬似乎是正确的，但没有 180° 环绕可以修复赤经，可能是因为它从春分开始。

我也曾尝试在 0,0 处使用观察器，但没有成功。

这可以使用 PyEphem、Skyfield 或 astropy 来完成吗？ PyEphem 中的人造卫星具有方便的亚纬度和亚长纬度属性，这似乎很奇怪，但对于天体来说却很难。

Answer 1

我终于明白了。有点。实际上我只是将 libastro 的相关部分移植到了 Python 中。请注意，此代码针对 Skyfield (be6c7296) 的当前 git 版本运行。

这里(gist version):

#!/usr/bin/env python3

from datetime import datetime, timezone
from math import atan, atan2, degrees, floor, pi, radians, sin, sqrt

from skyfield.api import earth, JulianDate, now, sun


def earth_latlon(x, y, z, time):
    """
    For an object at the given XYZ coordinates relative to the center of
    the Earth at the given datetime, returns the latitude and longitude
    as it would appear on a world map.

    Units for XYZ don't matter.
    """
    julian_date = JulianDate(utc=time).tt
    # see https://en.wikipedia.org/wiki/Julian_date#Variants
    # libastro calls this "mjd", but the "Modified Julian Date" is
    # something entirely different
    dublin_julian_date = julian_date - 2415020

    # the following block closely mirrors libastro, so don't blame me
    # if you have no clue what the variables mean or what the magic
    # numbers are because I don't either
    sidereal_solar = 1.0027379093
    sid_day = floor(dublin_julian_date)
    t = (sid_day - 0.5) / 36525
    sid_reference = (6.6460656 + (2400.051262 * t) + (0.00002581 * (t**2))) / 24
    sid_reference -= floor(sid_reference)
    lon = 2 * pi * ((dublin_julian_date - sid_day) *
                    sidereal_solar + sid_reference) - atan2(y, x)
    lon = lon % (2 * pi)
    lon -= pi
    lat = atan(z / sqrt(x**2 + y**2))

    return degrees(lat), degrees(-lon)


if __name__ == '__main__':
    print("2015-01-01 00:00:00:")
    time = datetime(2015, 1, 1, tzinfo=timezone.utc)
    x, y, z = earth(JulianDate(utc=time)).observe(sun).apparent().position.au
    print(earth_latlon(x, y, z, time))

    print("now:")
    time = datetime.now(timezone.utc)
    x, y, z = earth(JulianDate(utc=time)).observe(sun).apparent().position.au
    print(earth_latlon(x, y, z, time))

输出:

2015-01-01 00:00:00:
(-23.05923949080624, -179.2173856294249)
now:
(-8.384551051991025, -47.12917634395421)

如您所见，2015-01-01 00:00:00 的值与问题中的引用值匹配。不完全是，但对我来说已经足够了。据我所知，我的值(value)观可能会更好。

由于我对 libastro 代码中使用的未记录的魔数(Magic Number)一无所知，我无法使这项工作适用于地球以外的物体。

@BrandonRhodes:如果您有兴趣在 Skyfield 中使用此功能，请告诉我，然后我会尝试提出拉取请求。