PyEphem（パイ イーフェム）は、Pythonで天体の軌道や位置情報を扱うことができるライブラリです。内部的にはlibastroを利用しており、高精度な天体位置計算を手軽に行えます。
彗星（comet）や小惑星（asteroid）などの軌道計算のほか、SpaceXのStarlinkやOneWebなどの人工衛星の軌道を、TLE（Two-Line Element）データを用いて計算することができます。
本記事の内容
- PyEphemの基本的な使い方
- TLEデータを用いた人工衛星の位置計算
- 人工衛星の可視時刻の計算
- 地上軌跡の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
TLE（二行軌道要素形式）とは
TLE（Two-Line Element Set）は、人工衛星の軌道を記述するための標準的なフォーマットです。NORADが管理しており、CelesTrakから最新データを入手できます。
TLEは2行の数値データで構成され、以下の軌道要素を含みます。
$$ \text{TLE} \rightarrow \{i, \Omega, e, \omega, M, n\} $$
ここで、 * $i$: 軌道傾斜角 * $\Omega$: 昇交点赤経 * $e$: 離心率 * $\omega$: 近地点引数 * $M$: 平均近点角 * $n$: 平均運動（1日あたりの周回数）
PyEphemの基本
import ephem
import numpy as np
import matplotlib.pyplot as plt
# 観測地点の設定（東京・府中市）
observer = ephem.Observer()
observer.lat = '35.6694'
observer.lon = '139.4772'
observer.elevation = 60 # 標高 [m]
observer.date = '2026/2/19 12:00:00' # UTC
# 主要な天体の位置を計算
bodies = {
'Sun': ephem.Sun(),
'Moon': ephem.Moon(),
'Mars': ephem.Mars(),
'Jupiter': ephem.Jupiter(),
'Saturn': ephem.Saturn(),
}
print("=== 天体の位置（東京, 2026-02-19 12:00 UTC） ===")
for name, body in bodies.items():
body.compute(observer)
alt_deg = float(body.alt) * 180 / np.pi
az_deg = float(body.az) * 180 / np.pi
print(f"{name:10s}: 高度={alt_deg:+7.2f} deg, 方位角={az_deg:6.2f} deg")
人工衛星の可視時刻を計算する
PyEphemの next_pass 関数を使うと、地上の特定地点から人工衛星が可視になる時刻（AOS: Acquisition of Signal）、最大仰角時刻、不可視になる時刻（LOS: Loss of Signal）を計算できます。
import ephem
import numpy as np
# ISS の TLE データ
tle_name = "ISS (ZARYA)"
tle_line1 = "1 25544U 98067A 26049.50000000 .00016717 00000-0 10270-3 0 9003"
tle_line2 = "2 25544 51.6400 200.0000 0007000 90.0000 270.0000 15.50000000100001"
iss = ephem.readtle(tle_name, tle_line1, tle_line2)
# 観測地点（東京）
observer = ephem.Observer()
observer.lat = '35.6762'
observer.lon = '139.6503'
observer.elevation = 40
observer.date = '2026/2/19 00:00:00'
# 仰角の最低条件（度）
observer.horizon = '10'
# 次の5回分の可視パスを計算
print("=== ISS 可視パス（東京, 仰角10度以上） ===")
for i in range(5):
try:
info = observer.next_pass(iss)
rise_time, rise_az, max_time, max_alt, set_time, set_az = info
print(f"\nパス {i+1}:")
print(f" AOS (見え始め): {ephem.Date(rise_time)}, "
f"方位角={np.degrees(float(rise_az)):.1f} deg")
print(f" 最大仰角時刻: {ephem.Date(max_time)}, "
f"仰角={np.degrees(float(max_alt)):.1f} deg")
print(f" LOS (見え終わり): {ephem.Date(set_time)}, "
f"方位角={np.degrees(float(set_az)):.1f} deg")
# 次のパスを探すために時刻を進める
observer.date = set_time + ephem.minute
except ValueError:
print(f"パス {i+1}: 計算できませんでした")
observer.date += 1
複数衛星の地上軌跡を可視化する
import ephem
import numpy as np
import matplotlib.pyplot as plt
# 複数衛星の TLE データ（例）
satellites = {
"ISS": (
"1 25544U 98067A 26049.50000000 .00016717 00000-0 10270-3 0 9003",
"2 25544 51.6400 200.0000 0007000 90.0000 270.0000 15.50000000100001"
),
}
fig, ax = plt.subplots(figsize=(12, 6))
for name, (line1, line2) in satellites.items():
sat = ephem.readtle(name, line1, line2)
# 1軌道周期分の地上軌跡を計算
lats = []
lons = []
start = ephem.Date('2026/2/19 00:00:00')
# 約90分間（ISS 1周回分）
for i in range(900):
sat.compute(start + i * ephem.second * 6)
lats.append(np.degrees(float(sat.sublat)))
lons.append(np.degrees(float(sat.sublong)))
ax.scatter(lons, lats, s=1, label=name)
# 地図の枠
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel('Longitude [deg]')
ax.set_ylabel('Latitude [deg]')
ax.set_title('Satellite Ground Track')
ax.legend()
ax.grid(True, alpha=0.3)
# 赤道と主要緯線
ax.axhline(0, color='gray', linestyle='-', alpha=0.5)
for lat in [-66.5, -23.5, 23.5, 66.5]:
ax.axhline(lat, color='gray', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.show()
衛星の仰角プロファイル
特定のパスにおける衛星の仰角の時間変化をプロットしてみましょう。
$$ \text{El}(t) = \arcsin(\sin \varphi \sin \delta_s + \cos \varphi \cos \delta_s \cos h_s) $$
ここで、$\varphi$ は観測地点の緯度、$\delta_s$ は衛星の赤緯、$h_s$ は時角です。
import ephem
import numpy as np
import matplotlib.pyplot as plt
tle_name = "ISS (ZARYA)"
tle_line1 = "1 25544U 98067A 26049.50000000 .00016717 00000-0 10270-3 0 9003"
tle_line2 = "2 25544 51.6400 200.0000 0007000 90.0000 270.0000 15.50000000100001"
iss = ephem.readtle(tle_name, tle_line1, tle_line2)
observer = ephem.Observer()
observer.lat = '35.6762'
observer.lon = '139.6503'
observer.elevation = 40
observer.date = '2026/2/19 00:00:00'
observer.horizon = '0'
# 次のパスを取得
info = observer.next_pass(iss)
rise_time = info[0]
set_time = info[4]
# パス中の仰角・方位角を計算
duration = (set_time - rise_time) * 24 * 3600 # 秒
n_points = int(duration)
times = []
elevations = []
azimuths = []
for i in range(n_points):
t = rise_time + i * ephem.second
observer.date = t
iss.compute(observer)
times.append(i)
elevations.append(np.degrees(float(iss.alt)))
azimuths.append(np.degrees(float(iss.az)))
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 仰角の時間変化
axes[0].plot(times, elevations, 'b-', linewidth=2)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Elevation [deg]')
axes[0].set_title('ISS Elevation Profile')
axes[0].grid(True)
axes[0].axhline(10, color='r', linestyle='--', label='10 deg threshold')
axes[0].legend()
# 極座標でのスカイプロット
ax_polar = fig.add_subplot(122, projection='polar')
az_rad = np.radians(azimuths)
zenith = [90 - e for e in elevations]
ax_polar.plot(az_rad, zenith, 'b-', linewidth=2)
ax_polar.scatter(az_rad[0], zenith[0], c='green', s=100, zorder=5, label='AOS')
ax_polar.scatter(az_rad[-1], zenith[-1], c='red', s=100, zorder=5, label='LOS')
ax_polar.set_theta_zero_location('N')
ax_polar.set_theta_direction(-1)
ax_polar.set_ylim(0, 90)
ax_polar.set_title('Sky Plot')
ax_polar.legend(loc='lower right')
plt.tight_layout()
plt.show()
まとめ
本記事では、PyEphemを使ったPythonでの天体・人工衛星の軌道情報の扱い方を解説しました。
- PyEphemでは天体と人工衛星の位置を高精度に計算できる
- TLEデータを用いて人工衛星の位置・速度を任意の時刻で求められる
next_pass関数で衛星の可視パスを効率的に計算できる- 地上軌跡やスカイプロットの可視化により、衛星の運動を直感的に理解できる
次のステップとして、以下の記事も参考にしてください。