python - PyEphem:通过解析字符串表示形式转换负角时出现问题

Question

我一直在使用 PyEphem 来研究行星相对于地球上某个位置的运动。然而，我注意到有时从 PyEphem 得到的结果对于赤纬来说似乎是不连续的，而赤经似乎是连续的。这些图中每 30 分钟获取一次 RA 和 Dec 坐标。

我希望恒星体以连续的方式移动，但是当赤纬为负时，它看起来是不连续的。

有什么想法吗？如果这也有帮助的话，我可以发布我的脚本。

代码:

from ephem import *
import datetime
import re

def cSTI(number):
    degrees = int(number[0])
    minutes = int(number[1])
    seconds = float(number[2])
    full = degrees + (minutes/60) + (seconds/60/60)

    return round(full,2)

planets = [Sun(),Mercury(),Venus(),Moon(),Mars(),Jupiter(),Saturn(),Uranus(),Neptune(),Pluto()]
masses =  [1.989*10**30,.330*10**24,4.87*10**24,0.073*10**24,0.642*10**24,1898*10**24,568*10**24,86.8*10**24,102*10**24,0.0146*10**24]
earthMass = 5.97*10**24
gravitational_constant = 6.67408 * 10**-11
theYear = '2018'
theMonth = '9'
theDay = '8'
revere = Observer()
revere.lon = '-42.4084'
revere.lat = '71.0120'
currentDate = datetime.datetime(2018,11,30,12,0,0)
lowerDateLimit = datetime.datetime(2018,9,1,12,0,0)
revere.date = currentDate
print('DATE SUN(RA) SUN(DEC)    SUN(AZI)    SUN(GFORCE) MERCURY(RA) MERCURY(DEC)    MERCURY(AZI)    MERCURY(GFORCE) VENUS(RA)   VENUS(DEC)  VENUS(AZI)  VENUS(GFORCE)   MOON(RA)    MOON(DEC)   MOON(AZI)   MOON(GFORCE)    MARS(RA)    MARS(DEC)   MARS(AZI)   MARS(GFORCE)    JUPITER(RA) JUPITER(DEC)    JUPITER(AZI)    JUPITER(GFORCE) SATURN(RA)  SATURN(DEC) SATURN(AZI) SATURN(GFORCE)  URANUS(RA)  URANUS(DEC) URANUS(AZI) URANUS(GFORCE)  NEPTUNE(RA) NEPTUNE(DEC)    NEPTUNE(AZI)    NEPTUNE(GFORCE) PLUTO(RA)   PLUTO(DEC)  PLUTO(AZI)  PLUTO(GFORCE)   ')

while (currentDate> lowerDateLimit):
    print('%s   ' % (revere.date),end = ' ')
    planetIndex = 0;
    for planet in planets:
        planet.compute(revere)

        rightascension = str(planet.ra)
        declination = str(planet.dec)
        azimuth = str(planet.az)

        rightascension = re.split('[:]',rightascension)
        declination = re.split('[:]',declination)
        azimuth = re.split('[:]',azimuth )

        rightascension = cSTI(rightascension);
        declination = cSTI(declination);
        azimuth = cSTI(azimuth);
        GFORCE = gravitational_constant * ((masses[planetIndex]*earthMass)/(planet.earth_distance**2))
        print('%s   %s  %s  %s  ' % (rightascension,declination,azimuth,GFORCE),end = ' ')
        planetIndex+=1
    print()
    currentDate += datetime.timedelta(minutes=-30)
    revere.date = currentDate

Answer 1

我相信问题出在您从 ephem.Angle() (这是 ra、azi 等的对象)的手动转换中.) 到 float 。

<小时/>

编辑

特别是，出现问题的原因是在您的 cSTI() 函数中，当值为负数时，您应该减去(而不是添加)不同的值。

正确的实现如下所示:

import math

def cSTI(number):
    degrees = int(number[0])
    minutes = int(number[1])
    seconds = float(number[2])
    full = degrees + \
        math.copysign((minutes / 60), degrees) + \
        math.copysign((seconds / 60 / 60), degrees)
    return round(full, 2)

请注意，这是对代码进行的某种最小修改，以使其正常工作。我建议你写一些关于如何编写更好的代码的文档，从 PEP8 开始。另外，您应该避免在代码中使用魔数(Magic Number)，例如那些狂野的60。

如果我手动执行此操作，我还会避免使用不必要的正则表达式，并且实际上会从 ephem.Angle() 对象开始，例如:

import math

MIN_IN_DEG = SEC_IN_MIN = 60

def ephem_angle_to_float(angle):
    degrees, minutes, seconds = [float(s) for s in str(angle).split(':')]
    value = abs(degrees) + \
        minutes / MIN_IN_DEG + \
        seconds / SEC_IN_MIN / MIN_IN_DEG
    return math.copysign(value, degrees)

并且可以在代码中直接调用 planet.ra 或 planet.dec。

但我不会手动执行此操作。见下文。

<小时/>

但好消息是，您无需自己手动计算，您只需将其转换为 float 即可获得以弧度为单位的值，如 official documentation 中所示。 .

以下代码是您的改进版本，它以列表的形式生成表格数据(可以使用 Python 标准库中的 csv 轻松将其转换为 CSV)。进行了一些修改:

• 天体和质量现在更加明确
• ephem 对象是动态创建的
• G 常数现在取自 scipy.constants(从最新的 CODATA 中获取)
• 标签和数据是动态创建的
• 尽可能避免显式索引
• 开始日期和结束日期有些颠倒，因为这与问题无关

请注意，此处未执行任何转换，因为这会使我们丢失信息。

这是代码:

import numpy as np
import scipy as sp

import ephem
import scipy.constants
import datetime

celestial_masses = {
    'sun': 1.989e30,
    'mercury': 0.330e24,
    'venus': 4.870e24,
    'earth': 5.970e24,
    'moon': 0.073e24,
    'mars': 0.642e24,
    'jupiter': 1898e24,
    'saturn': 568e24,
    'uranus': 86.8e24,
    'neptune': 102e24,
    'pluto': 0.0146e24, }
celestials = {
    name: (eval('ephem.{}()'.format(name.title())), mass)
    for name, mass in celestial_masses.items() if name != 'earth'}
gg = sp.constants.gravitational_constant

observer = ephem.Observer()
# Revere, Massachusetts, USA
observer.lon = '-42.4084'
observer.lat = '71.0120'

start_datetime = datetime.datetime(2018, 9, 1, 12, 0, 0)
end_datetime = datetime.datetime(2018, 11, 30, 12, 0, 0)

values = 'ra', 'dec', 'azi', 'gforce'
labels = ('date',) + tuple(
    '{}({})'.format(name, value)
    for name in celestials.keys()
    for value in values)
data = []
observer.date = start_datetime
delta_datetime = datetime.timedelta(minutes=30)
while (observer.date.datetime() < end_datetime):
    row = [observer.date]
    for name, (body, mass) in celestials.items():
        body.compute(observer)
        row.extend([
            body.ra, body.dec, body.az,
            gg * ((mass * celestial_masses['earth']) / (body.earth_distance ** 2))])
    data.append(row)
    observer.date = observer.date.datetime() + delta_datetime

要转换为 CSV(使用 ra、dec、az 和 gforce 的浮点值)，可以这样做:

import csv

filepath = 'celestial_bodies_revere_MA.csv'
with open(filepath, 'w') as file_obj:
    csv_obj = csv.writer(file_obj)
    csv_obj.writerow(labels)
    for row in data:
        # : use default string conversion
        # csv_obj.writerow(row)

        # : a possible conversion to float for all but the `date`
        csv_obj.writerow([
            float(x) if i != labels.index('date') else x
            for i, x in enumerate(row)])

此外，这里是一些绘图代码，显示负值的问题已经消失:

import matplotlib as mpl
import matplotlib.pyplot as plt

x = [row[labels.index('date')].datetime() for row in data]
fig, axs = plt.subplots(4, 1, figsize=(16, 26))
for i, value in enumerate(values):
    pos = i
    for name in celestials.keys():
        y = [row[labels.index('{}({})'.format(name, value))] for row in data]
        axs[pos].plot(x, y, label=name)
        axs[pos].set_title(value)
        axs[pos].legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout()

产生: