StarLinkの軌道を地図に描写する

StarLinkのネタが続きますが、今回はStarLinkの衛星の位置と軌道を地図にプロットしてみたいと思います。

軌道情報を入手する

軌道情報(TLE)は衛星の軌道を決定する力学的パラメータです。この情報があれば衛星の軌道を決定し、衛星の位置を計算することができます。 StarLink衛星のTLEは以下のサイトから入手可能です。

https://celestrak.org/NORAD/elements/gp.php?GROUP=starlink&FORMAT=tle

skyfieldで地図にプロットする

今回はpythonを使用します。天文関連の計算を気軽に行えるskyfieldモジュールを使用します。

from skyfield.api import Topos, load, EarthSatellite,wgs84
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import math
import numpy as np
from numpy import linalg as LA

tle_fille = 'starlink_tle.txt'  #ダウンロードしたTLEファイルのパス

with open(tle_fille) as f:
    l_strip = [s.rstrip() for s in f.readlines()]

ts = load.timescale()
t0 = ts.utc(2023, 7, 11, 6, 20, 40)

lats=[]
lons=[]

for i in range(0, len(l_strip), 3):
    starlink= EarthSatellite(l_strip[i+1], l_strip[i+2], l_strip[i], ts)
    geocentric = starlink.at(t0)  # 地球中心を基準としたx, y, zの位置を取得
    subpoint = geocentric.subpoint()  # 緯度、経度、高度情報に変換
    
    lats.append(subpoint.latitude.degrees)
    lons.append(subpoint.longitude.degrees)


fig=plt.figure()
ax=fig.add_axes([0.1,0.1,5.0,5.0])
# setup mercator map projection.
m = Basemap(llcrnrlon=-180.,llcrnrlat=-85.,urcrnrlon=180.,urcrnrlat=85.,\
            resolution='l',projection='merc')

m.scatter(m(lons, lats)[0],m(lons, lats)[1], marker = 'o', color='r', zorder=10)

m.drawcoastlines()
m.fillcontinents()

plt.show()

プロットしてみた図が上です。ほとんどの衛星が北緯/南緯55度あたりで折り返しています。また列になっている衛星はStarLinkトレインと呼ばれ、夜間であれば連なっているのが肉眼で見ることができます。 時間が経つにつれて徐々に衛星の間隔が開き、連なっているのは見れなくなります。

参考文献

qiita.com