（ほかのアニメーション画像も記事の下のほうにあります。Safariで画像がカクカクする場合は画像を開いてご覧ください。）
はじめに
私たちはGPS衛星（をはじめとした、ガリレオ、みちびき、などのGNSS衛星）との位置関係を用いて自分の居場所を正確に知ることができます。ということは、それらの人工衛星自体が今どこを飛んでいるのかもきっと取得できるはずです。
少し調べてみると、実はGPS/GNSSの衛星というのは、衛星の現在位置そのものを送信しているわけではなく、軌道を表すためのパラメータ（＝軌道要素）を送信していることがわかりました。利用者側は、その軌道要素をもとに衛星の現在位置を手元で計算するわけですね。
そして、人工衛星の軌道というのは当然ながら様々な機関に監視されているため、軌道要素のデータは、GPSデバイスなどで直接受信せずとも、CelesTrak などのサイトから誰でもダウンロードできることも知りました（※ただし、実際の測位で使うものではなく、精度も全く違うことに注意が必要です）。
……というところまで分かったところで、「軌道のパラメータが得られるなら、簡易なアニメーションを自分でさくっと描けるはず！」と思いついて、実際に試してみることにしました。
軌道のデータを読む
CelesTrakの NORAD GP Element Sets Current Data のページに行き、例えば「GNSS」の所をクリックすると このようなテキストデータ がダウンロードできます。
GPS BIIR-2 (PRN 13)
1 24876U 97035A 23041.68353357 .00000027 00000+0 00000+0 0 9998
2 24876 55.5468 146.3758 0065839 52.4873 308.1231 2.00564062187432
(...以下同じような形式のデータが続く...)
ファイルを眺めると、GPSをはじめ、EUのガリレオ (GSAT)、中国の北斗 (BEIDOU)、日本のみちびき (QZS)など、様々な国や地域の測位用衛星が含まれています。なお、軌道のようすは日々変っていくため、このデータも随時更新されていることに注意してください。
このデータは TLE（2行軌道要素形式）と呼ばれる形式になっています。PythonでTLEデータを扱えるライブラリを探したところ、Skyfield というライブラリが扱いやすそうだったので、今回はこれを採用することにしました。早速 TLE を読み込んでみます。
>>> import skyfield.api
# TLEを読み込む
>>> sats = skyfield.api.load.tle_file("https://celestrak.org/NORAD/elements/gnss.txt", reload=True)
# 衛星が何個含まれているか？
>>> len(sats)
156
# 0番目の衛星をみてみる
>>> sats[0]
<EarthSatellite GPS BIIR-2 (PRN 13) catalog #24876 epoch 2023-01-26 17:25:39 UTC>
156個あるそれぞれの衛星が EarthSatellite というオブジェクトで表現されています。
衛星の位置を取得する
ここから衛星の位置を得るのは簡単で、EarthSatellite.at(t) というメソッドを呼ぶだけで指定した時刻の衛星の位置を計算してくれます（裏側ではSGP4という軌道計算アルゴリズムが走るようです）。
>>> ts = skyfield.api.load.timescale()
>>> sat = sats[0]
# 衛星の現在の位置を取得してみる
>>> sat.at(ts.now())
<Geocentric ICRS position and velocity at date t center=399 target=-124876>
# 座標を取得してみる
>>> sat.at(ts.now()).position
<Distance [-1.48650828e-04 2.47026029e-05 9.16709615e-05] au>
# メートル単位で欲しい
>>> sat.at(ts.now()).position.m
array([-22229143.44338772, 3673091.13254689, 13733781.9073863 ])
# 速度も得られる
>>> sat.velocity.m_per_s
array([ 1611.64563397, -2842.84923151, 2128.79653024])
あっさりと3次元座標らしきものが手に入りました。
軌道要素は日々更新されています。計算対象の時点が軌道要素データの更新時点から離れるほど精度が落ちていくことに注意が必要です。
ここで位置として返される Geocentric オブジェクトは、地球の重心を原点としつつ地球の自転と一緒に回転しない座標系 (GCRS; Geocentric Celestial Reference System) での位置を表しています（なお、軸は 国際天文基準座標系 (ICRS) と同じ方向を向いているようです）。自転している地球を宇宙から眺めているイメージですね。
ここまでで、「指定した時刻の、地球を中心とした座標」が簡単に手に入ることがわかったので、もうアニメーションが描けてしまいそうです。
3次元でプロットしてみる
得られた座標をそのままPythonの matplotlib を使って3次元のアニメーションプロットにしてみます。（matplotlibでアニメーションを描く方法は今回の範囲を超えてしまうので省略しますが、興味のある方は FunctionAnimation などで調べてみてください）。
図の描き方として：
- matplotlib の3Dプロット (
projection="3d") を使う。 - GNSSのシステムごとに色を塗りわける。
- 衛星の位置を点で描く (scatter plot)。おまけとして速度ベクトルも描く (quiver plot)。
- 軌跡を描くため、5分ごとの位置を線でつないで描く。
というようにしてみたところ、下図のような結果になりました。なんだかそれらしい軌道が描けました！
繰り返しますが、これは地球自体に固定されていない座標系 (GCRS) で見ています。赤道のあたりをまわっているものは静止衛星ですが、地球の回転に固定されていない座標系なので、位置は固定されずに地球の自転と同じ速さで（同期して）周回しています。
地球の回転と固定してみる
上の例は地球に固定されていない座標系でしたが、地球に固定された座標系からも見てみましょう。そのような座標を Skyfield で得るには以下のようにします。ちなみにここで使っている座標系は ITRS/ITRF（国際地球基準座標系）というものです。このような座標系を ECEF (Earth-Centered, Earth-Fixed) などと呼んだりもするようです。
from skyfield.framelib import itrs
# ITRSでの座標を得る
xyz = sat.at(t).frame_xyz(itrs)
# 速度も得たい場合は
(xyz, velocity) = pos.frame_xyz_and_velocity(itrs)
コードをほんの少し変えるだけで、地球に固定された座標が得られました。この座標系でも同じようにアニメーションを描いてみましょう。
わかりやすくするために日本の「みちびき (QZSS)」を太く描いてみました。地球の回転に固定された座標系で見るとちゃんと南北非対称の8字形で日本のあたりを周回しています（個人的にはちょっと感動しました）。ちなみにその隣にいる赤い8字形軌道の衛星は中国の北斗です。
地図も付けたい
3次元プロットだけだと少しさみしい気もしたので、地図上にもプロットしてみました。matplotlib で地図のプロットをする場合、いまは cartopy というライブラリが定番のようです。
衛星の位置を地図上にプロットするには、人工衛星の座標 (x, y, z) を、相当する地理経緯度 (lon, lat) に変換してあげる必要があります。Skyfield で経緯度 (WGS84) を得るには以下のようにします。
from skyfield.api import wgs84
pos = sat.at(t)
(lat, lng) = wgs84.latlon_of(pos) # 衛星の位置に相当する経緯度を取得
こうして得られた経緯度を使い、catropy の cartopy.crs.PlateCarree() 投影法で地図を付けたものが、記事冒頭にも載せたアニメーションです。このように航跡を地上の位置で描いたものを ground trace や ground track などと呼ぶそうです。
ちなみに、cartopy の Nightshade という機能を使って夜の場所に影をつけています。
おまけ： 某衛星コンステレーションを描いてみる
CelesTrak などのサイトでは、GNSSだけでなく、宇宙ステーションや、気象衛星など様々な人工衛星の軌道要素が公開されています。
最後に、最近話題の某衛星コンステレーションのデータを読み込んでみました。
