🛰️

Juliaで人工衛星の軌道をプロットする

2021/03/16に公開

Julia

はじめに

前回はSatelliteToolboxを使用して特定の地点からいつ見えるかを計算しました。
今回は反対に人工衛星が地球上をどのように動いているか地図上にプロットしてみたいと思います。

実行環境は前回と同様にJupyterLab上でJuliaを動かしています。
(v1.5.3)

SatelliteToolboxの導入、TLEの取得までは前回の記事を参考にしてください。

TLEの取得

SatelliteToolboxの導入

ISSの軌道を計算してみる

ISSのTLEは前回と同様のものを利用します。
軌道プロパゲーターも同じくSGP4です。

軌道取得関数

前回は地上局の可視時間計算用の便利な関数がありました。
軌道計算も同じく便利な関数が用意されているのでそれを利用することにしましょう。

ground_trace(orbp::OrbitPropagator{N}, eop_data::Union{Nothing, EOPData_IAU1980, EOPData_IAU2000A} = nothing; ECI = TEME(), ECEF = PEF(), span = 1.0) where N

最低限必要なパラメーターは前回も使用した軌道プロパゲーターの1つだけです。

eopについてはあまり良くわかっていませんがなくても問題ないようです。
ECI、ECEFも前回同様ですがデフォルト値が入っているため省略しても良いです。

span

一つだけ設定したいオプションがspanです。
これは何周分のデータを取得したいか設定する値です。
デフォルトとして1が入っていますが今回は3周分取得してみましょう。

戻り値

戻り値として(lat, lon)で一組のタプルの配列が[(lat1, lon1), (lat2, lon2), ...]のようなデータが取得できます。

計算

ということで計算は単純に下記のように記述すればよいです。

trace = ground_trace(orbp, span=3.0)

1673-element Array{Tuple{Float64,Float64},1}:
 (0.8831537534586672, 2.5300798162991742)
 (0.8805436426046043, 2.546561971564309)
 (0.8777887803463595, 2.562933057547235)
 (0.8748909133714676, 2.57918855135821)
 (0.8718518495589704, 2.595324190695959)
 (0.8686734530902965, 2.6113359852551055)
 (0.8653576395639468, 2.627220227016527)
 (0.8619063711391004, 2.642973472759355)
 (0.8583216517316619, 2.658592551862235)
 (0.8546055222844771, 2.6740745730411657)
 (0.8507600561316158, 2.689416903615516)
 (0.8467873544746825, 2.7046171742462484)
 (0.8426895419871605, 2.719673282813654)
 ⋮
 (0.9020625582535058, 1.1032062661225848)
 (0.9011714635834893, 1.1205498771684652)
 (0.9001216080990879, 1.1378493506420506)
 (0.8989137469305198, 1.1550977345622893)
 (0.8975487428688709, 1.172288222381454)
 (0.8960275634015791, 1.1894141793898383)
 (0.8943512774311715, 1.2064691416476594)
 (0.8925210517027451, 1.2234468400459435)
 (0.8905381469674017, 1.2403412232402733)
 (0.888403913910091, 1.257146452670606)
 (0.8861197888713672, 1.2738569227979326)
 (0.88368728939308, 1.2904672798007306)

プロット

一秒ごとの軌道のリストが得られたのでプロットしてみてみます。
プロットするときに緯度経度を反対にしたほうが楽なのと、ラジアンから度数にしたほうがわかりやすいので変換しておきます。

trace = map(x->(rad2deg(x[2]), rad2deg(x[1]), trace)

この状態で一度プロットしてみましょう。

plot(trace, label=nothing, xlims=(-180, 180), ylims=(-90, 90))

するとこのような感じに出力されると思います。
(取得したTLEの情報により軌道の位置は異なります)
これはISSの軌道は左から右に移動しているのですが、経度が180から-180に戻るときにそのままプロットされてしまっています。
この現象についてかんたんに対応できる方法ないかなーと思ってましたが、Julia関連の情報をよく出してくださっている黒木玄さんよりTwitterで情報をいただきました。

実直にデータを分割するような情報を与えてあげるしかなさそうです。
上記投稿のスレッドでいくつかやり方に関してやり取りしてますが、今回は端っこに行く度にデータを分ける方法で考えてみます。

データ分割

ということで次のような関数を定義します。

function split_data(input; th=100)
    output = [similar(input, 0)]
    tmp = similar(input, 0)
    for i in 2:length(input)
        x, y = input[i].-input[i-1]
        if abs(x)>th || abs(y)>th
            push!(output, tmp)
            tmp = similar(input, 0)
        end
        push!(tmp, input[i])
    end
    push!(output, tmp)
    output
end

引数としてinputに先程生成したtraceを渡してあげれば良いようにしています。
thはしきい値で隣同士の緯度経度でどのくらい値が飛んだら分割したと認識するかを決めてます。
とりあえず人工衛星が1秒で100度も移動することはないのでこのくらいで十分だと思います。（何なら10くらいでもいいですが）

さてこれで分割した値を改めてプロットしてみます。
複数の配列に分割されているのでfor文でプロットを繰り返します。

trace_split = split_data(trace)
plt = plot() #プロットの初期化
for i in 1:length(trace_split)
    plot!(trace_split[i], label=nothing, linecolor =:blue)
end
plt

背景に地図を表示する

ここまでできたら軌道に関しては問題ないと思いますが、背景がなにもないのは味気ないので地図を表示してみましょう。

地図データを用意する前に地図と地図投影法に関して軽く説明します。
まず地球は丸いので、球体上の位置を示す地形情報はそのまま2次元平面に持ってくることはできません。
何らかの投影を行い2次元にマッピングしなければならないのですが、その方法にはいくつか種類があり目的により使い分けたりします。
例えば国同士の面積を比較したい時は、面積が一定に保たれるが緯度の線が樽状に歪んだ地図(正積図法)を使用します(ハンメル図法等)が、移動距離等を見たい時は距離を等しくしたいので正距図法を使用したりします。
詳しくはWikipediaやその他資料を調べてみてください。

正距円筒図法

今回は緯度経度の情報を簡単に扱える形式が良いですよね。
そのような場合は正距円筒図法を使うのがかんたんです。
名前の通り球体を円筒と見立てて投影したものです。
緯度・経度がそれぞれ縦軸、横軸の長さと対応しているので変換処理が必要なくそのまま描画することができます。
その代わり円筒に見立てているので北極・南極側が引き伸ばされており面積としては不正確なものになっています(ロシアがとても大きく見える)
その点に目をつぶれば楽ですし、衛星軌道の表示方法としてもよく使用されているので問題はないでしょう。

地図データ

地図データを用意しないといけませんが自由に使えて手軽な地図としては、Wikipediaの正距円筒図法のページに載ってある地図かなと思います。

Equirectangular-projection.jpg

このページの情報を見るとパブリックドメインとなっているので自由に使用して問題なさそうですのでこの画像をダウンロードしてJuliaのプロジェクトから読み込める場所に保存してください。

地図データ表示

画像データをJuliaで使用するにはImageパッケージを導入する必要があります。
またJPEGを扱うには追加でImageMagicも必要みたいなので一緒に宣言しておきます。
画像読み込みはload(filename)を使用すればいいので次のように記述します。

Pkg.add("Images")
Pkg.add("ImageMagick")
using Images

img = load("Equirectangular-projection.jpg")

画像を表示するにはimgをplotに渡してあげれば良いのですが、そのままだと画像データの場合自動で縦軸の上が0で値が大きくなるについれて下に行きます。
画像データとしては自然で扱いやすいのですが、今回は緯度経度として扱いたいのでちょっと都合が悪いです。
なのでデータを渡すときに img[end:-1:1, :]のようにして反転させ、更にplotの引数としてyflip = falseを指定して正しい方向に表示できるようにします。
地図を緯度経度フルで表示できるようにxlims, ylimsも指定して次のように記述すると表示できます。

plt = plot([-180, 180], [-90, 90], img[end:-1:1, :], yflip = false, size=(900, 450), xlims=(-180, 180), ylims=(-90, 90))

これで全部揃ったので改めて処理全体を載せます。

trace = map(x->(rad2deg(x[2]), rad2deg(x[1])), ground_trace(orbp, span=3.0))
split_trace = split_data(trace)
img = load("Equirectangular-projection.jpg")
plt = plot([-180, 180], [-90, 90], img[end:-1:1, :], yflip = false, size=(900, 450), xlims=(-180, 180), ylims=(-90, 90))
for i in 1:length(tmp)
    plot!(split_trace[i], label=nothing, linecolor=:red)
end
plt

得られる画像は次のようになります。

きれいに描画できていますね！

おわりに

長々と説明してしまいましたが、処理本体はほとんど何もすることなくとてもシンプルに記述できますね。
前回のground_station_access同様自由度は低いですが、かんたんに使える関数が用意されていて少しずつ学ぶには便利でいいですね。
今後も引き続き使い方を調べていこうかなと思います。
次はもっと高度な計算を行う部分を書いていければと思います。

何か質問や指摘がございましたらご気軽にコメントください。