Starlinkの衛星軌道を地図上にプロットしてみる

# はじめに

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

# 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
```

[こちら](https://www.celestrak.com/NORAD/elements/gp.php?GROUP=active&FORMAT=tle)から適当にコピペしています。
衛星の名称は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　です。

https://rhodesmill.org/pyephem/index.html

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

https://rhodesmill.org/pyephem/tutorial.html

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

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

ようです。

## プロットしてみる

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

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

## 結果

![out.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294040/343e5907-daee-5af0-da24-c3f891e06fcc.png)

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

https://satellitemap.space/?constellation=starlink

https://www.reddit.com/r/Starlink/comments/jl28r4/starlink_visualization_satellitemapspace/?utm_source=share&utm_medium=web2x&context=3

# まとめ

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


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

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

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

https://xephem.github.io/XEphem/Site/xephem.html

https://www.clearskyinstitute.com/xephem/help/xephem.html

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

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

https://qiita.com/Ushinji/items/76d4901658bce767df36

(2022/04/03)

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