はじめに

なんだか衛星軌道をプロットしてみたくなり、してみました。衛星軌道情報は公開されています。衛星軌道情報を読み込み、時刻を与えると衛星位置を計算するライブラリもあります。世界地図上に任意の緯度経度の点をプロットすることもできます。ということで、衛星軌道をプロットできます。

PyephemでTLEを読んで軌道をプロット

衛星軌道情報

これはTLE(two line elements)というものが有名です。

こちらでSpace X のもので例えば、以下のものがありました。

STARLINK-61             
1 44249U 19029Q   22092.64473324  .00156203  00000+0  44400-2 0  9992
2 44249  52.9765  64.4520 0001512 266.1363  93.9465 15.35088908157866

こちらから適当にコピペしています。
衛星の名称はSTARLINK-61です。

行番号１の行を読みます。

要素名	値
人工衛星番号	44249
分類	U（分類無し）
国際識別番号	19(打ち上げ年の下二桁)
国際識別番号	029(打ち上げ年における打ち上げ数)
国際識別番号	Q
エポック年	22(西暦の下2桁)
エポック	092.64473324(年通算日+その日における00時からの経過時間)

軌道情報は行番号２の行です。ケプラーの軌道六要素というらしいですが、軌道長半径，離心率，軌道傾斜角，近地点引数，昇交点赤経，近地点通過時刻から構成されるそうです。これらをまとめて紹介するのは別として、ここでは先のTLEからとってきます。TLEの読み方はこちらなどに情報があります。

https://gportal.jaxa.jp/gpr/assets/mng_upload/GCOM-W/TLE.pdf

要素名	値
人工衛星番号	44249
軌道傾斜角（度）	52.9765
昇交点赤径（度）	64.4520
離心率	0.0001512
近地点引数（度）	266.1363
平均近点離角（度）	93.9465
平均運動（周回数/日）	15.35088908

Pyephem

衛星軌道情報を読んで、日時を与えると衛星位置を計算してくれる便利なツールはPyehem　です。

ここのtutorial に書いてある通りにすれば計算できます。

軌道パラメータはXEPHEMのformat というもので読み込めるようですが、ここでは先程紹介したTLEを用いました。

ephem.EarthSatellite インスタンスをreaddle()でTLEを読んで作成する
このインスタンスにcompute()で時刻を指定すると、その時刻の位置をシミュレートする

ようです。

プロットしてみる

以下のコードでプロットした図を作成できました。コードに出てくるyear と doy(day of year)は先に紹介したＴＬＥの行番号１の情報から取得しています。

plot_orbit_starlink.py

import ephem
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from datetime import datetime, timedelta
 
tle_start_link="""STARLINK-61             
1 44249U 19029Q   22092.64473324  .00156203  00000+0  44400-2 0  9992
2 44249  52.9765  64.4520 0001512 266.1363  93.9465 15.35088908157866"""

year, doy = 2022, 92.64473324

dt0 = datetime(year, 1, 1) + timedelta(days=doy)
time_list = [(dt0 + timedelta(hours=i/60)) for i in range(0,60*12)]

m = Basemap()
m.drawcoastlines()

tle0, tle1, tle2 = tle_start_link.split("\n")
satellite = ephem.readtle(tle0, tle1, tle2)
for dt in time_list:
    item = dt.strftime("%Y/%m/%d %H:%M:%S")
    satellite.compute(f'{item}')
    latitude = satellite.sublat / ephem.degree
    longitude = satellite.sublong / ephem.degree
    x1,y1 = m(longitude, latitude)
    m.plot(x1, y1, "m.", markersize=1)
 
plt.savefig('out.png')
#plt.show()

結果

プロットしてみて、この衛星は北極や南極まで行かないので「本当に合っているの？」と心配になりましたが、下記をみるとそれで良さそうです。それでもいささか不安もありますが。^^;

まとめ

TLEから軌道情報を読み、衛星軌道を世界地図上にプロットできるようになりました。もっと使い勝手を良くしたいなぁ。と思いますが、とりあえず、いろいろ楽しめそうです。

参考にさせていただきたページ

みちびきを軌道の図示、参考にさせて頂きました。
https://korintje.com/archives/77

今回は関係がありませんでしたが、ここら辺のXEPHEMは気になります。

TLEとか自分でparseしてみるかな。
http://celestrak.com/NORAD/documentation/tle-fmt.php
地図ももっときれいなのを使いたい。

こちらでは pyorbital というライブラリを使っていた。Cesium の開発についても触ってみたいな。

(2022/04/03)

参考情報でpyorbital のページを追加 (2022/04/11)