使用GMT绘制星下点图

写在前面

这是2023年第一篇文章，今天要介绍的是使用GMT绘制星下点图，在卫星导航实验重，要求我们绘制北斗卫星的星下点图，我们需要知道卫星的轨迹，就需要计算出卫星的坐标，并进行相应的坐标转换，获得了坐标之后，需要使用GMT进行星下点图的绘制。

一、获得卫星坐标

可以在我的其他的博客中找到计算卫星坐标的代码，此处放置链接，请自行取用。

【python】读取卫星星历（RENIX 3.04）进行卫星位置的计算（北斗卫星专题）

其大致原理是，通过读取卫星星历相应的卫星轨道信息，对相应的数据块进行读取并计算出在历元时刻的卫星位置。

二、对卫星坐标进行坐标转换

由于计算获得的卫星坐标是以高斯空间直角坐标系（CGCS2000或其他坐标系），我们需要转化为BLH坐标，才能够使用GMT来进行卫星的星下点图。以下是转换坐标的相关类的代码。

1.建立椭球类，采用CGCS2000坐标系

class Earth():
    # 建立椭球类，，采用CGCS2000中国大地坐标系数据
    def __init__(self):  # 给定椭球数据
        # 长半轴
        self.a = 6378137.0
        # 短半轴
        self.b = 6356752.31414
        # 扁率
        self.f = 1 / 298.257222101
        # 第一偏心率的平方
        self.e1s = m.pow(0.0818191910428, 2)
        # 第二偏心率的平方
        self.e2s = self.e1s / (1 - self.e1s)

2.建立转换类（包括XYZ与BLH的互相转换方法）

class Transform():
    # 建立转换类
    def __init__(self):  # 导入椭球数据
        self.earth = Earth()

    def BLH_XYZ(self, B, L, H):  # 大地坐标转换为空间直角坐标
        # 将角度转换为弧度
        self.b = m.radians(B)
        self.l = m.radians(L)
        # 计算辅助函数
        self.W = m.sqrt(1 - self.earth.e1s * m.pow(m.sin(self.b), 2))
        # 转换为空间直角坐标
        self.N = self.earth.a / self.W
        self.X = (self.N + H) * m.cos(self.b) * m.cos(self.l)