Python in Astronomy (4.数据的处理与可视化（一））

Astro_ChaoXu

已于 2024-02-23 11:00:58 修改

阅读量994

点赞数 25

文章标签： python 开发语言

于 2024-02-22 21:06:14 首次发布

本文链接：https://blog.csdn.net/m0_51889082/article/details/136081540

版权

NumPy是使用Python进行科学计算的基础包，支持大量的维度数组与矩阵运算。我们将结合NumPy和Astropy库，计算天体在天球上的位置。此外我们还可引入Matplotlib库从数组数据中生成可视化的图像。

2.1球面天文学

在天文学中，恒星的位置坐标可以由两个角坐标来确定，而恒星与地球的径向距离是第三个坐标，但是并不是所有的天体都需要知道这个距离，在实际运用中，把所有天体沿镜像投射到一个以地球为中心的球体上，这个球体也就是所谓的天球。

对于地球上的观测者来说，天体的位置坐标取决于地理纬度和经度，并随一天中时间的变化而变化，这两个角坐标也被称为高度角和方位角。如果分别在北京和东京同时观测同一天体，不难发现，天体会对应不同的高度角和方位角。为了得到独立于观察者的位置坐标，可使用相对于固定参考方向的角坐标。天文学中最常用的是赤道坐标系，在赤道坐标系中，一个参考方向是由地球的自转轴给出的，从恒星到天赤道平面的角距离称为赤纬，用δ表示。另一个参考方向是由赤道平面和地球绕太阳轨道运动平面（黄道面）的交点所定义。如图所示，黄道与赤道有两个交点，每年春季太阳所在的那一点称为春分点，另一点是秋分点。赤道坐标系中的第二个坐标称为赤经，用α表示。赤经从春分点开始沿着天赤道向东度量。

通过上图可以更清楚地明白什么是赤纬？什么是赤经？约定俗成，赤纬从赤道向天北极和天南极度量，正负各90度，向北为正向南为负。赤经从春分点开始，沿着与天体东升西落相反的方向度量，以时角、时分、时秒的单位代替度、分、秒的单位，从 $0^h - 24^h$ 。换算关系为：

$1 ^h = 15^{\circ}, \quad 1^m = 15{}', \quad 1^s = 15{}''$

2.2太阳赤纬

在赤道坐标系中，恒星的赤纬是恒定的，但是太阳在赤道坐标系中的位置在一年的时间里会发生变化。这是地球自转轴相对于垂直于黄道的方向倾斜的结果，这个角度被称为黄赤夹角 $\large{\boldsymbol{\varepsilon_0}}$ ，约等于23.44度，太阳赤纬在一年中的变化量，可以用下面近似公式来描述：

$\delta _ {\odot} = \displaystyle -arc \sin \left[ \sin \large{\boldsymbol{\varepsilon_0}} \cos \left( \frac {360^{\circ} }{365.24} (N+10 \right ) \right]$

其中N为从1月1日开始的天数之差。所以一年的第一天对应N = 0，最后一天对应N = 364（除非是闰年）。分式360/365.24度是地球角坐标的日变化量（假设是一个圆形轨道），这也是地球每天绕太阳轨道运动的角速度ω。太阳在春分或秋分（赤道和黄道的交点）赤纬为0，在夏至日或冬至日赤纬达到 $\pm \large{\boldsymbol{\varepsilon_0}}$ ，在那里地球的旋转轴朝向太阳或远离太阳。确切的日期每年都有所不同。例如，在2020年的3月20日和9月22日，6月20日和12月21日。

下面的Python代码根据这个公式计算给定的N值的赤纬 $\delta _ {\odot}$ ：

import math

N = 171 #冬至日
omega = 2*math.pi/365.24 # 角速度
ecl = math.radians(23.44) # 黄赤交角

# 太阳赤纬的近似表达式
delta = -math.asin(math.sin(ecl)*math.cos(omega*(N+10)))
print("declination = {:.2f} deg".format(math.degrees(delta)))

结果接近于23.44度的预期值。在上面代码中我们使用了math库中的正弦、余弦和反正弦函数，这些函数的参数必须用弧度来指定，角速度由表达式omega = 2*math.pi/365.24 给出，黄赤交角则通过 math.radians()函数转化为弧度制。

为了计算第二个至日和对应的赤纬，我们需要计算第8行的代码中对应的N值。我们可以重新给N赋值，然后再次运行代码得到计算结果，虽然这对于一些不同的值肯定是可行的，但对于许多值来说，它将变得过于乏味。理想情况下，我们想计算出几次太阳的赤纬，这可以通过使用被称为数组的数据结构来实现。数组是具有相同类型的数据元素的有序集合。下面是一个例子：

import numpy as np

# 2020年的春分，夏至，秋分，冬至
N = np.array([79, 171, 265, 355])

与其他编程语言相比，数组不是Python原生的。它们是在numpy模块中定义的。我们这里使用的是别名np导入numpy库，函数np.array（）接受括号中的值列表，将这些值作为元素创建一个数组。与数组一样，列表也是数据元素的有序集合。它们的区别与联系如下所示：

数组操作的语法与列表相似
数组 (Numpy array) 要求元素的数据类型被预设且一致，列表 (List) 无此要求
数组的存储是一段连续的内存空间，列表不是

由np.array（）返回的数组被分配给变量N，它的主要属性可以通过以下打印语句来显示：

print(N)
print(N.size)
print(N.dtype)

由此我们可以看到，N有四个元素（元素的数量是通过N.size属性获得的），四个元素分别是整数79、171、256和355。数据类型可以使用.dtype属性来显示。

我们如何处理数组中的值？例如，N[0]是指数组的第一个元素，也就是数字79。通过方括号中的一个整数可以引用数组中的单个元素，这个数指定了该元素在数组中的位置，被称为索引。在Python中，第一个元素的索引为0，第二个元素的索引为1，以此类推。您还可以从数组的末尾开始索引元素，最后一个元素的索引为−1，最后二个元素的索引为−2，以此类推。

print(N[1])
print(N[-3])

通过以上两行代码均可输出数字171，您可以继续改变索引，并观察能得到什么结果，这有助于帮助您更好的理解数组的索引。总的来说，数组中的每个元素都是相同的数据类型，且可以通过其索引来标识。

将N定义为一个数组后，我们可以通过复制上面的代码中的delta表达式，并用数组N替换表达式中的N来计算太阳的赤纬。

delta = -np.arcsin(math.sin(ecl) * np.cos(omega*(N+10)))
print(np.degrees(delta))

上面代码可对数组中的每一个元素进行计算，并将结果储存在一个名为delta的新数组中，值得注意的是前面还有一个delta，但是会给后面这个表达式所覆盖，所以输出的是覆盖之后的内容，其4个元素分别对应着它们的赤纬值。这个结果并不完美（等分处的赤纬应该是零），但相当接近，我们毕竟只使用了一个近似的公式。

为了更好地了解上面的代码是如何工作的，我们可将其扩展为几个步骤，并使用一个名为tmp的临时数组来处理计算的中间结果：

tmp = N+10
print(tmp)
print(tmp.dtype)

tmp = omega*tmp
print(tmp)
print(tmp.dtype)

tmp = math.sin(ecl)*np.cos(tmp)
print(tmp)

delta = -np.arcsin(tmp)
print(np.degrees(delta))

如您所见，我们可以对数组执行算术运算，可以是单个数组和单个数字运算，也可以是两个具有相同形状的数组运算。例如，在上述代码的第一部分中，其中的数字10会分别和数组N中的每一个元素相乘，然后分配给数组tmp。

下一步是与ω（地球每天绕太阳轨道运动的角速度）的乘法运算，数组tmp中的每一个元素都会和ω相乘，然后再次分配给数组tmp。由于ω的值是一个浮点数，所以tmp数组会自动从数据类型整数转换为浮点数：

NumPy函数np.cos（）可用于计算tmp数组中每个值的余弦值，而math.sin（）只能取一个单值参数，即ecl（黄道的倾角）。运算得到的值将会分配给数组tmp。得到的结果如下，

用np.arcsin计算每个元素的反正弦，最终得到赤纬值：

由上，我们可以总结以下的numpy基本操作：

现在，假设我们想要很好地格式化元素，就像我们在本节开始时在引入数组之前所做的那样。通过索引该元素，可以实现数组的特定元素的格式化打印。例如，

print("declination = {:.2f} deg".format(math.degrees(delta[1])))

要想打印所有的元素，我们可以使用循环遍历该数组，下面的示例显示了一个可用的循环：

for val in delta: 
    print("declination = {:6.2f} deg".format(math.degrees(val)))

该循环从数组中的第一个元素开始，将其值赋给变量val，然后执行循环体。上面的示例中的循环只包含一个打印语句。执行此语句后，循环继续使用下一个元素，以此类推，直到达到数组的结束（最高索引），得到的结果如下所示：

然而，这仍不是最好的方案。我们虽然打印输出了每一个赤纬值，但是并不知道它所对应的时间，为了打印天数和相应的赤纬值，我们需要同时遍历N和delta的元素。一种解决方案是利用Python的enumerate() 函数定义一个计数器：

print("i day delta [deg]")
for i, val in enumerate(delta):
    print("{1:d} {2:3d} {0:8.2f}".format(math.degrees(val), i, N[i]))

表达式中的占位符{1：d}表示第二个参数（冒号前的索引1）将作为整数插入到此位置，下一个占位符{2：3d}引用第三个参数（冒号前面的索引2），最后一个占位符引用第一个参数。如果没有指定，则一个萝卜一个坑，按顺序填入。

那么是否可以编写一个循环，同时迭代两个数组，而不使用计数器i呢？实际上，这就是zip()函数的功能，它从两个或多个相同长度的数组中聚合具有相同索引的元素。在这种情况下循环变量是一个包含每个数组中的一个元素的元组。当然，这里的元组又是一种新的数据类型，但是其使用频率并不是很多，这里不详细介绍。目前所需要知道的是，以下代码示例中的元组行包含N的一个元素和相应的delta元素，形成了我们要打印的表的一行。元组元素像数组元素一样被索引，row[0]指的是当天，row[1]指的是当天的赤纬值。

print("day delta [deg]")
for row in zip(N,delta):
    print("{0:3d} {1:8.2f}".format(row[0],math.degrees(row[1])))

除了索引列，我们得到与上述相同的输出：

尝试将代码的第三行改成print(row)，看看会得到什么？

print("day delta [deg]")
for row in zip(N,delta):
    print(row)

方括号用于列表和数组，而元组则用方括号括起来。最重要的区别是，元组是不可变的，也就是说，不可能添加、更改或删除元素。虽然可以通过设置N[3] = 364来修改数组N，但为元组中的单个元素赋值是不可以的，您只能用一个新的元组来覆盖整个元组。值得注意的是上面输出的结果并没有错误，只是用弧度制表示的，在格式化输出结果前，请不要忘了考虑输出物理量的单位。

通过以上示例，我们可自行的更改数组N中的天数，比如下面就在数组N中添加了6个日期，并且得到了其赤纬值。

参考文献

[1]. W. Schmidt and M. Völschow, Numerical Python in Astronomy and Astrophysics , Undergraduate Lecture Notes in Physics,
[2]. 苏宜编著. 2019. 天文学新概论（第五版）[M]. 北京：科学出版社