matplotlib Basemapは、Python経由で2次元の地図を扱うためのライブラリです。
主に、海洋学者や気象学者などが利用することを想定しており、様々な座標参照系（Coordinate Reference System, CRS）を相互変換する機能や、陸地や国境線などを描写する機能を備えています。
本記事の内容
- Basemapでできることの概要
- 基本となるBasemapオブジェクトと投影法
- 地図上へのデータプロット
- 衛星の地上軌跡を描画する実践例
Basemapのインストール
Basemapはcondaまたはpipでインストールできます。
conda install -c conda-forge basemap
# または
pip install basemap
なお、BasemapはCartopyへの移行が推奨されていますが、シンプルな地図描画では依然として便利なツールです。ここではCartopyを使った代替実装も併記します。
座標参照系（CRS）の基礎
地球上の位置を平面に投影するには、座標参照系が必要です。地球は球体（正確には回転楕円体）であるため、平面への投影には必ず歪みが生じます。
主要な投影法
| 投影法 | 特徴 | 用途 |
|---|---|---|
| 正距円筒図法（Plate Carree） | 緯度・経度をそのまま直交座標に | 全球データの概観 |
| メルカトル図法 | 角度を正確に保つ | 航海用海図 |
| ランベルト正積方位図法 | 面積を正確に保つ | 統計地図 |
| 正距方位図法 | 中心からの距離が正確 | 航空路 |
| 極射影 | 極域の観測に最適 | 極域研究 |
各投影法には固有の歪みがあり、用途に応じて使い分けます。これは球面を平面に写像する以上、角度・面積・距離を同時に保つことは不可能であるためです（ガウスの驚異の定理）。
Cartopyを使った地図描画
Basemapの後継として推奨されるCartopyを使った実装方法を紹介します。
import numpy as np
import matplotlib.pyplot as plt
try:
import cartopy.crs as ccrs
import cartopy.feature as cfeature
USE_CARTOPY = True
except ImportError:
USE_CARTOPY = False
if USE_CARTOPY:
# Cartopyを使った全球地図
fig, ax = plt.subplots(figsize=(12, 6),
subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_global()
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.BORDERS, linewidth=0.3, linestyle='--')
ax.gridlines(draw_labels=True, linewidth=0.3, alpha=0.5)
ax.set_title('World Map (Plate Carree Projection)')
plt.tight_layout()
plt.show()
else:
# Cartopyが使えない場合のシンプルな代替
fig, ax = plt.subplots(figsize=(12, 6))
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel('Longitude [deg]')
ax.set_ylabel('Latitude [deg]')
ax.set_title('World Map (Simple)')
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)
plt.tight_layout()
plt.show()
投影法の数学
正距円筒図法（Plate Carree）
最もシンプルな投影法で、緯度 $\varphi$ と経度 $\lambda$ をそのまま直交座標に対応させます。
$$ x = \lambda, \quad y = \varphi $$
メルカトル図法
正角図法として知られ、緯度 $\varphi$ を次のように変換します。
$$ x = \lambda, \quad y = \ln \tan\left(\frac{\pi}{4} + \frac{\varphi}{2}\right) $$
この変換により、地図上の任意の点で角度が保存されます。ただし、高緯度に行くほど面積が拡大されるため、面積の比較には不向きです。
球面上の2点間距離（ハヴァサイン公式）
球面上の2点 $(\varphi_1, \lambda_1)$ と $(\varphi_2, \lambda_2)$ の距離は、ハヴァサイン公式で求められます。
$$ d = 2R \arcsin\sqrt{\sin^2\left(\frac{\varphi_2 – \varphi_1}{2}\right) + \cos \varphi_1 \cos \varphi_2 \sin^2\left(\frac{\lambda_2 – \lambda_1}{2}\right)} $$
import numpy as np
def haversine(lat1, lon1, lat2, lon2, R=6371.0):
"""ハヴァサイン公式による2点間距離 [km]"""
lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
dlat = lat2 - lat1
dlon = lon2 - lon1
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
return 2 * R * np.arcsin(np.sqrt(a))
# 東京 -> ニューヨーク
d = haversine(35.6762, 139.6503, 40.7128, -74.0060)
print(f"東京 -> ニューヨーク: {d:.0f} km")
# 東京 -> ロンドン
d = haversine(35.6762, 139.6503, 51.5074, -0.1278)
print(f"東京 -> ロンドン: {d:.0f} km")
衛星の地上軌跡を地図に描画する
import numpy as np
import matplotlib.pyplot as plt
# 衛星軌道パラメータ
Re = 6371.0 # 地球半径 [km]
h = 500.0 # 軌道高度 [km]
inc = np.radians(97.4) # 軌道傾斜角（太陽同期軌道）
mu = 398600.4418 # 地球重力定数 [km^3/s^2]
a = Re + h
T = 2 * np.pi * np.sqrt(a**3 / mu) # 軌道周期 [s]
# 地上軌跡の計算（3周回分）
omega_e = 2 * np.pi / 86164.1 # 地球の自転角速度 [rad/s]
omega_s = 2 * np.pi / T # 衛星の角速度 [rad/s]
t = np.linspace(0, 3 * T, 10000)
lat = np.degrees(np.arcsin(np.sin(inc) * np.sin(omega_s * t)))
lon_inertial = np.degrees(np.arctan2(
np.sin(omega_s * t) * np.cos(inc),
np.cos(omega_s * t)
))
lon = np.mod(lon_inertial - np.degrees(omega_e * t) + 180, 360) - 180
# 経度の不連続点を処理
lon_diff = np.abs(np.diff(lon))
lon[:-1][lon_diff > 180] = np.nan
fig, ax = plt.subplots(figsize=(14, 7))
ax.scatter(lon, lat, c=t/60, cmap='viridis', s=1)
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel('Longitude [deg]')
ax.set_ylabel('Latitude [deg]')
ax.set_title(f'Ground Track (h={h:.0f}km, inc={np.degrees(inc):.1f}deg, 3 orbits)')
ax.grid(True, alpha=0.3)
# 主要な緯度線
for lat_line in [-66.5, -23.5, 0, 23.5, 66.5]:
ax.axhline(lat_line, color='gray', linestyle='--', alpha=0.3)
cb = plt.colorbar(ax.collections[0], ax=ax)
cb.set_label('Time [min]')
plt.tight_layout()
plt.show()
print(f"軌道周期: {T/60:.1f} min")
print(f"1日の周回数: {86400/T:.1f}")
まとめ
本記事では、matplotlib BasemapおよびCartopyを使った地図情報の扱い方を解説しました。
- Basemap/Cartopyは座標参照系の変換と地図描画を簡単に行えるライブラリである
- 投影法は用途に応じて正距円筒・メルカトル・正距方位などから選択する
- ハヴァサイン公式で球面上の2点間距離を正確に計算できる
- 人工衛星の地上軌跡の可視化に活用できる
次のステップとして、以下の記事も参考にしてください。