【Python】Skyfieldで衛星軌道予測

# はじめに
以前、[こちら](https://qiita.com/grinpeaceman/items/23af7733962268bd900c)の記事で、[Pyorbital](https://pyorbital.readthedocs.io/en/feature-moon-phase/)を使った衛星の軌道計算についてご紹介しました。
しかし、Pyorbitalのライセンスは[GPL](https://ja.wikipedia.org/wiki/GNU_General_Public_License)なので、GPLのパッケージを使うのはちょっと・・・
という方向けに、今回は[Skyfield](https://pypi.org/project/skyfield/)を使って衛星の軌道計算をやってみました。

# Skyfieldとは
Pythonの天文学計算パッケージです。
地球の周りの軌道にある星、惑星、衛星の位置などを計算することができます。
ライセンスは[MITライセンス](https://ja.wikipedia.org/wiki/MIT_License)です。

# 環境
+ Windows 10
+ Python 3.7.3
+ skyfield 1.23

# インストール
`$ pip install skyfield`

# 使ってみた
## [TLE](https://ja.wikipedia.org/wiki/2%E8%A1%8C%E8%BB%8C%E9%81%93%E8%A6%81%E7%B4%A0%E5%BD%A2%E5%BC%8F#:~:text=2%E8%A1%8C%E8%BB%8C%E9%81%93%E8%A6%81%E7%B4%A0%E5%BD%A2%E5%BC%8F%EF%BC%88%E3%81%AB%E3%81%8E%E3%82%87%E3%81%86%E3%81%8D%E3%81%A9%E3%81%86%E3%82%88%E3%81%86%E3%81%9D,%E5%BD%A2%E5%BC%8F%E3%81%AE%E3%83%95%E3%82%A9%E3%83%BC%E3%83%9E%E3%83%83%E3%83%88%E3%81%A7%E3%81%82%E3%82%8B%E3%80%82)読み込み

衛星のTLEを用意。今回は[CELESTRAK](https://celestrak.com/)から地球観測衛星[Aqua](https://ja.wikipedia.org/wiki/Aqua_(%E4%BA%BA%E5%B7%A5%E8%A1%9B%E6%98%9F))の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

```

```python: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

## 衛星の位置情報取得

```python: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


## パス時刻の取得

```python: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に設定しとかないと怒られます。