はじめに
以前、こちらの記事で、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
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
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
衛星の位置情報取得
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
パス時刻の取得
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に設定しとかないと怒られます。