はじめに
なんだか衛星軌道をプロットしてみたくなり、してみました。衛星軌道情報は公開されています。衛星軌道情報を読み込み、時刻を与えると衛星位置を計算するライブラリもあります。世界地図上に任意の緯度経度の点をプロットすることもできます。ということで、衛星軌道をプロットできます。
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です。
行番号1の行を読みます。
要素名 | 値 |
---|---|
人工衛星番号 | 44249 |
分類 | U(分類無し) |
国際識別番号 | 19(打ち上げ年の下二桁) |
国際識別番号 | 029(打ち上げ年における打ち上げ数) |
国際識別番号 | Q |
エポック年 | 22(西暦の下2桁) |
エポック | 092.64473324(年通算日+その日における00時からの経過時間) |
軌道情報は行番号2の行です。ケプラーの軌道六要素というらしいですが、軌道長半径,離心率,軌道傾斜角,近地点引数,昇交点赤経,近地点通過時刻から構成されるそうです。これらをまとめて紹介するのは別として、ここでは先のTLEからとってきます。TLEの読み方はこちらなどに情報があります。
要素名 | 値 |
---|---|
人工衛星番号 | 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)は先に紹介したTLEの行番号1の情報から取得しています。
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)