Pythonでの測位衛星スカイプロット作成
はじめに
測位衛星からの電波受信のみで自らの座標を知るためには、少なくとも4衛星からの測位信号が必要です。利用者位置において、緯度、軽度、楕円体高、時刻の4変数が未知なためです。
測位にはより多くの衛星が観測できることが好ましいです。このような天空の測位衛星の位置を表示するPythonコードを作成しました。
スカイプロット
例えば、内閣府が公開している衛星配置表示アプリGNSS Viewを用いると、任意地点・任意時刻での測位衛星の見える位置がわかります。これは、2023年11月23日 23:20:00 UTC（日本時間の24日 08:20:00）に、広島にて観測できるであろう欧州Galileo衛星の配置を、GNSS Viewにてプロットしたスクリーンショットの一部です。
この図面は、自らの位置において、測位衛星の見える方角と仰角を表すもので、スカイプロット（skyplot）と呼ばれます。
千葉工業大学の鈴木太郎先生が公開されているオープンソースのGNSS-Radarは、インターネット接続なしにブラウザ上に測位衛星の位置を表示できます。私も測位衛星の位置を知るためにGNSS-Radarを利用しています。
このように、測位衛星の位置を知る手段はあります。
一方、レポート作成のためには、PDF形式のグレイスケール画像が好ましく、また、特定衛星をハイライト表示したり、観測条件を図面上に表示するためにも、自らの手でスカイプロットを作成したいです。また、LaTeXなどにこのような図面を読み込ませるためには、PNGファイルやJPGファイルのようなビットマップファイルではなく、ベクトル形式にて画像収録したPDFファイルが好ましいです。
自らのコードでスカイプロットを作成するため、衛星軌道を入手します。
2行軌道要素
衛星測位を行うためには、その衛星の正確な軌道を表すエフェメリス（ephemeris）が必要です。一方、スカイプロット表示のためには大まかな衛星位置がわかればよく、このような用途には、しばしば2行軌道要素（TLE: two line element）が用いられます。これは、衛星軌道表現に必要な軌道傾斜角などを2行のテキスト形式にて表したものです。TLEの詳細は、きどうようそのひみつにわかりやすく記載されています。
このようなTLEはCelesTrakなどから入手できます。厳密には、衛星軌道は刻々と変化するものなのかもしれませんが、軌道修正中でなければ、TLEは大きくは変化しないと考えられます。エフェメリスの寿命は2時間ですが、TLEは、一度ダウンロードすれば、しばらく利用できます。例えば、先日打ち上げに成功した先進レーダ衛星 だいち4号のTLEは次のように表されます。
ALOS-4 (DAICHI-4)
1 60182U 24123A 24229.54940804 .00000500 00000+0 73855-4 0 9998
2 60182 97.9215 325.3457 0001395 98.1122 28.2230 14.79473854 6861
ここでは、欧州のGalileo全衛星のTLEをCelesTrakから入手してみます。
このLatest Data欄にある書類マークをクリックすると、ブラウザ上にTLEが表示されますので、テキストファイル上にコピーペーストします。すべてのGalileo衛星に対して、この作業をくり返します。
Pythonによるスカイプロット作成
次に、TLEをもとにして、任意座標点からの衛星の見える方角や仰角を計算します。座標計算や表示のために、Pythonの非標準モジュールskyfield、numpy、matplotlibを利用しました。
pip install skyfield numpy matplotlib
ここで作成したコードは、辞書satellites1に格納された利用可能衛星を表す丸印で表示する一方、satellites2に格納された利用不可衛星を表すバツ印で表示します。さらに、利用できない衛星番号はグレーアウトされるようにします。ここでは、利用可能衛星として、E03、E07、E08、E13、E15を挙げています。EはGalileo衛星を表す欧州（Europe）のEです。TLEの記述順は任意です。観測点は広島（座標は北緯34.23、東経132.27）で、日時は2023年11月22日 23:30:00 UTCです。
skyplot.py
#! /usr/bin/env python
# -*- coding: utf-8 -*-
FILENAME = 'skyplot'
FONTSIZE = 16
# GAL TLE can be obtained from:
# https://celestrak.org/norad/elements/table.php?GROUP=galileo&FORMAT=tle
from skyfield.api import Topos, load, wgs84
from skyfield.sgp4lib import EarthSatellite
import numpy as np
from matplotlib import pyplot as plt
import datetime
# TLE data
satellites1 = {
"E03": ["1 41860U 16069B 24075.03714384 .00000000 00000+0 00000+0 0 9997",
"2 41860 55.0417 123.9214 0003428 309.7138 50.3288 1.70474893 45607"],
"E07": ["1 41859U 16069A 24075.76931632 -.00000003 00000+0 00000+0 0 9996",
"2 41859 55.0449 123.9045 0005163 290.7537 69.2607 1.70474715 45362"],
"E08": ["1 41175U 15079B 24077.45505889 -.00000001 00000+0 00000+0 0 9994",
"2 41175 55.3797 123.8317 0004223 301.5343 58.4922 1.70474689 51340"],
"E13": ["1 43567U 18060D 24073.22747821 -.00000102 00000+0 00000+0 0 9991",
"2 43567 57.3051 4.0560 0004572 309.8424 50.1177 1.70475162 35092"],
"E15": ["1 43564U 18060A 24073.15310948 -.00000102 00000+0 00000+0 0 9996",
"2 43564 57.3028 4.0542 0005044 303.0742 56.8777 1.70475060 35093"],
}
satellites2 = {
"E11": ["1 37846U 11060A 24070.59038482 -.00000094 00000+0 00000+0 0 9990",
"2 37846 57.1269 4.2009 0004048 0.9151 359.0813 1.70475681 76996"],
"E12": ["1 37847U 11060B 24075.36853573 -.00000108 00000+0 00000+0 0 9993",
"2 37847 57.1277 4.0716 0005128 340.5417 117.7805 1.70475226 77083"],
"E19": ["1 38857U 12055A 24077.08212499 -.00000003 00000+0 00000+0 0 9991",
"2 38857 55.3621 124.1355 0005270 267.6966 92.3107 1.70473587 70991"],
"E20": ["1 38858U 12055B 24075.29159274 -.00000001 00000+0 00000+0 0 9992",
"2 38858 55.3626 124.1864 0002996 263.1608 96.8769 1.70473746 70965"],
"E18": ["1 40128U 14050A 24076.38068143 -.00000063 00000+0 00000+0 0 9995",
"2 40128 49.7505 306.7902 1607231 147.6212 223.3476 1.85519539 63054"],
"E14": ["1 40129U 14050B 24077.17831918 -.00000057 00000+0 00000+0 0 9995",
"2 40129 49.7668 305.8335 1605771 148.4561 222.2616 1.85520304 65236"],
"E26": ["1 40544U 15017A 24073.30254099 -.00000102 00000+0 00000+0 0 9991",
"2 40544 56.9222 4.0745 0004144 291.6220 68.3342 1.70475448 55129"],
"E22": ["1 40545U 15017B 24074.14320105 -.00000105 00000+0 00000+0 0 9990",
"2 40545 56.9283 4.0665 0003375 281.1733 78.7899 1.70475327 18189"],
"E24": ["1 40889U 15045A 24076.84456596 .00000028 00000+0 00000+0 0 9998",
"2 40889 55.2656 244.3185 0005887 44.3383 315.6720 1.70473651 52993"],
"E30": ["1 40890U 15045B 24075.30543699 .00000030 00000+0 00000+0 0 9994",
"2 40890 55.2642 244.3563 0003840 46.2777 313.7129 1.70473690 52990"],
"E09": ["1 41174U 15079A 24077.23534912 -.00000002 00000+0 00000+0 0 9999",
"2 41174 55.3811 123.8397 0004688 303.6832 56.3397 1.70474783 50768"],
"E02": ["1 41549U 16030A 24076.99092123 .00000028 00000+0 00000+0 0 9995",
"2 41549 55.4130 244.3111 0003383 18.6658 341.3102 1.70473716 48653"],
"E01": ["1 41550U 16030B 24076.10899448 .00000030 00000+0 00000+0 0 9998",
"2 41550 55.4111 244.3364 0001211 357.6320 2.3289 1.70473751 48631"],
"E04": ["1 41861U 16069C 24077.75112543 -.00000001 00000+0 00000+0 0 9990",
"2 41861 55.0457 123.8527 0005491 274.6415 85.3637 1.70473447 45518"],
"E05": ["1 41862U 16069D 24074.96434658 .00000001 00000+0 00000+0 0 9994",
"2 41862 55.0430 123.9236 0004700 275.5406 84.4789 1.70474973 45591"],
"E21": ["1 43055U 17079A 24076.04030420 .00000030 00000+0 00000+0 0 9994",
"2 43055 55.3883 244.1605 0000736 22.1381 337.8261 1.70474387 38952"],
"E25": ["1 43056U 17079B 24076.91865947 .00000028 00000+0 00000+0 0 9998",
"2 43056 55.3883 244.1362 0001720 316.3645 43.5851 1.70474520 38986"],
"E27": ["1 43057U 17079C 24073.03300014 .00000027 00000+0 00000+0 0 9997",
"2 43057 55.3924 244.2419 0002160 350.7256 9.2233 1.70474360 38903"],
"E31": ["1 43058U 17079D 24076.77253818 .00000029 00000+0 00000+0 0 9993",
"2 43058 55.3864 244.1366 0001974 300.9948 58.9487 1.70474451 38996"],
"E33": ["1 43565U 18060B 24075.13123991 -.00000108 00000+0 00000+0 0 9993",
"2 43565 57.3054 4.0028 0003358 282.1792 77.7846 1.70475082 35161"],
"E36": ["1 43566U 18060C 24073.59351374 -.00000103 00000+0 00000+0 0 9995",
"2 43566 57.3052 4.0464 0004542 314.0609 45.9024 1.70475319 35115"],
"E34": ["1 49809U 21116A 24074.42529437 -.00000106 00000+0 00000+0 0 9990",
"2 49809 57.2961 3.9185 0001376 281.4869 183.9787 1.70475239 14122"],
"E10": ["1 49810U 21116B 24075.39474642 -.00000108 00000+0 00000+0 0 9998",
"2 49810 57.2968 3.8907 0001282 308.3744 144.5609 1.70475412 14162"],
}
# Load timescale and observer position
ts = load.timescale()
t = ts.utc(2023, 11, 22, 23, 30, 0)
observer = Topos(latitude_degrees=34.23, longitude_degrees=132.27)
# Initialize satellite objects
satellite_objects1 = {sat: EarthSatellite(satellites1[sat][0], satellites1[sat][1], sat, ts) for sat in satellites1}
satellite_objects2 = {sat: EarthSatellite(satellites2[sat][0], satellites2[sat][1], sat, ts) for sat in satellites2}
# Plotting
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'},
figsize=(6, 6))
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.set_ylim(0, 90) # Altitude 0° (horizon) to 90° (zenith)
for sat in satellite_objects1:
difference = satellite_objects1[sat] - observer
topocentric = difference.at(t)
alt, az, distance = topocentric.altaz()
if alt.degrees > 0: # Only plot if above the horizon
r = 90 - alt.degrees
theta = np.radians(az.degrees)
ax.plot(theta, r, 'o', c='black')
ax.text(theta, r, f'{sat}', fontsize=FONTSIZE, ha='right', va='bottom')
for sat in satellite_objects2:
difference = satellite_objects2[sat] - observer
topocentric = difference.at(t)
alt, az, distance = topocentric.altaz()
if alt.degrees > 0: # Only plot if above the horizon
r = 90 - alt.degrees
theta = np.radians(az.degrees)
ax.plot(theta, r, 'x', c='gray')
ax.text(theta, r, f'{sat}', fontsize=FONTSIZE, ha='right', va='bottom', c='gray')
ax.set_rmax(90)
ax.set_rticks([0, 30, 60, 90]) # Altitude rings
ax.set_rlabel_position(-22.5) # Move radial labels
ax.set_rgrids([0, 30, 60, 90], labels=["90", "60", "30", "0"], fontsize=FONTSIZE)
ax.grid(True)
ax.set_thetalim([0, 2*np.pi])
ax.set_thetagrids(np.rad2deg(np.linspace(0, 2*np.pi, 5)[1:]),
labels=["E", "S", "W", "N"], fontsize=FONTSIZE)
plt.savefig(FILENAME + '.pdf')
#plt.savefig(FILENAME + '.svg')
#plt.savefig(FILENAME + '.tiff', dpi=600)
#plt.savefig(FILENAME + '.jpg', dpi=600)
#plt.savefig(FILENAME + '.png', transparent=True, dpi=600)
plt.show()
./skyplot.pyにてこのコードを実行すると、画面上にスカイプロットが表示され、そのPDFファイルも出力されます。ソースコードの修正により、出力形式をSVG、TIFF、JPG、または、PNG形式にしたり、観測点や観測時刻を連続変化させたり、様々な情報の追加表示させることも可能です。
まとめ
CelesTrakから2行軌道要素TLEを入手して、自らのPythonコードでスカイプロットを作成しました。手元にこのようなコードがあると、任意地点での、任意衛星のスカイプロットを作成することができ、様々な検証のためにも便利です。