はじめに

以前、こちらの記事で、Pyorbitalを使った衛星の軌道計算についてご紹介しました。
しかし、PyorbitalのライセンスはGPLなので、GPLのパッケージを使うのはちょっと・・・
という方向けに、今回はSkyfieldを使って衛星の軌道計算をやってみました。

Skyfieldとは

Pythonの天文学計算パッケージです。
地球の周りの軌道にある星、惑星、衛星の位置などを計算することができます。
ライセンスはMITライセンスです。

環境

Windows 10
Python 3.7.3
skyfield 1.23

インストール

$ pip install skyfield

使ってみた

TLE読み込み

衛星のTLEを用意。今回はCELESTRAKから地球観測衛星AquaのTLEを入手。

aqua.txt

AQUA                    
1 27424U 02022A   20203.40846355  .00000031  00000-0  17013-4 0  9997
2 27424  98.2101 143.5537 0000452 148.6720 321.9438 14.57110132968816

tle_read.py

from skyfield.api import load


satellite = load.tle_file("./aqua.txt")[0]  # TLE読み込み
sgp4_model = satellite.model

print("衛星名: ", satellite.name)
print("元期: ", satellite.epoch.utc_jpl())
print("衛星番号: ", sgp4_model.satnum)
print("エポック年: ", sgp4_model.epochyr)
print("エポック日: ", sgp4_model.epochdays)
print("エポックのユリウス日: ", sgp4_model.jdsatepoch)
print("ndot: ", sgp4_model.ndot)
print("nddot: ", sgp4_model.nddot)
print("弾道抗力係数 B*: ", sgp4_model.bstar)
print("軌道傾斜角[rad]: ", sgp4_model.inclo)
print("昇交点赤経[rad]: ", sgp4_model.nodeo)
print("離心率: ", sgp4_model.ecco)
print("近地点引数[rad]: ", sgp4_model.argpo)
print("平均近点角[rad]: ", sgp4_model.mo)
print("平均運動[rad/min]: ", sgp4_model.no_kozai)

実行結果

衛星名: AQUA
元期: A.D. 2020-Jul-21 09:48:11.2507 UT
衛星番号: 27424
エポック年: 20
エポック日: 203.40846355
エポックのユリウス日: 2459051.5
ndot: 9.39326507149726e-13
nddot: 0.0
弾道抗力係数 B*: 1.7013e-05
軌道傾斜角[rad]: 1.714089603712883
昇交点赤経[rad]: 2.505484718420184
離心率: 4.52e-05
近地点引数[rad]: 2.5948159055250097
平均近点角[rad]: 5.618979316382121
平均運動[rad/min]: 0.06357842341892297

衛星の位置情報取得

satellite_position.py

from skyfield.api import load


satellite = load.tle_file("./aqua.txt")[0]  # TLE読み込み

ts = load.timescale()  # タイムスケールを生成
now = ts.now()
geocentric = satellite.at(now)  # 地球中心を基準としたx, y, zの位置を取得
subpoint = geocentric.subpoint()  # 緯度、経度、高度情報に変換

print("時刻: ", now.utc_jpl())
print("緯度[degree]", subpoint.latitude.degrees)
print("経度[degree]", subpoint.longitude.degrees)
print("高度[m]", subpoint.elevation.m)

実行結果

時刻: A.D. 2020-Jul-22 00:49:04.0105 UT
緯度[degree]: 29.340257418205187
経度[degree]: 16.29265032283678
高度[m]: 705898.890496908

パス時刻の取得

get_pass_time.py

from datetime import datetime, timedelta, timezone

from skyfield.api import load, Topos


satellite = load.tle_file("./aqua.txt")[0]  # TLE読み込み

ts = load.timescale()  # タイムスケールを生成

# 観測地点(東京タワー)の位置情報を設定
tokyo_tower_pos = Topos(
    latitude_degrees=35.66, longitude_degrees=139.75, elevation_m=333
)

start = datetime.now(timezone.utc)
end = start + timedelta(days=1)

# 24時間以内に観測点(東京タワー)から衛星が見える時間を計算(仰角5度以上)
# パス開始時刻、最大仰角時刻、パス終了時刻の情報が取得できる
t, events = satellite.find_events(
    tokyo_tower_pos, ts.utc(start), ts.utc(end), altitude_degrees=5.0
)

difference = satellite - tokyo_tower_pos
for ti, event in zip(t, events):
    name = ("パス開始: ", "最大仰角: ", "パス終了: ")[event]
    print(name, ti.utc_strftime("[UTC] %Y %b %d %H:%M:%S"))
    if event == 1:
        topocentric = difference.at(ti)
        alt, _, _ = topocentric.altaz()  # 観測点から観た衛星の仰角を取得
        print("最大仰角[degree]: ", alt.degrees)
    if event == 2:
        print()

実行結果

パス開始:  [UTC] 2020 Jul 22 03:30:00
最大仰角:  [UTC] 2020 Jul 22 03:35:40
最大仰角[degree]:  58.23825328597088
パス終了:  [UTC] 2020 Jul 22 03:41:22

パス開始:  [UTC] 2020 Jul 22 05:09:57
最大仰角:  [UTC] 2020 Jul 22 05:13:44
最大仰角[degree]:  11.91191755398313
パス終了:  [UTC] 2020 Jul 22 05:17:32

パス開始:  [UTC] 2020 Jul 22 15:33:05
最大仰角:  [UTC] 2020 Jul 22 15:37:48
最大仰角[degree]:  18.970133046864028
パス終了:  [UTC] 2020 Jul 22 15:42:30

パス開始:  [UTC] 2020 Jul 22 17:10:10
最大仰角:  [UTC] 2020 Jul 22 17:15:37
最大仰角[degree]:  37.15713012973206
パス終了:  [UTC] 2020 Jul 22 17:21:03

まとめ

Pyorbitalでできることは、Skyfieldでも大体できそう。
注意点としては、計算に使うタイムスケールオブジェクトをdatetimeから生成するときに、タイムゾーンをUTCに設定しとかないと怒られます。