map_flight.py

import folium
import numpy as np

#指定した２点間を球体である地球を考慮し最短経路で結ぶスクリプト（日付変更線を横切る）
# 球面三角法を使用して大円を計算する関数
def great_circle_points(point1, point2, num_points):
    lat1, lon1 = np.radians(point1)
    lat2, lon2 = np.radians(point2)

    # 球面三角法による角度（ラジアン）を計算
    angles = np.linspace(0, 1, num_points)
    distance = np.arccos(np.sin(lat1) * np.sin(lat2) + np.cos(lat1) * np.cos(lat2) * np.cos(lon2 - lon1)) # 球面三角法の余弦定理の式変形（２点と北極点の3点を結んだ球面三角形を考える）

    # 大円上の各点の緯度経度を計算
    points = []
    for angle in angles:
        x = np.sin((1 - angle) * distance) / np.sin(distance) * np.cos(lat1) * np.cos(lon1) + np.sin(angle * distance) / np.sin(distance) * np.cos(lat2) * np.cos(lon2)
        y = np.sin((1 - angle) * distance) / np.sin(distance) * np.cos(lat1) * np.sin(lon1) + np.sin(angle * distance) / np.sin(distance) * np.cos(lat2) * np.sin(lon2)
        z = np.sin((1 - angle) * distance) / np.sin(distance) * np.sin(lat1) + np.sin(angle * distance) / np.sin(distance) * np.sin(lat2)
        lat = np.arctan2(z, np.sqrt(x ** 2 + y ** 2))
        lon = np.arctan2(y, x)

        # 経路を180度の子午線を横切るように修正
        if lon - lon1 < -np.pi:
            lon += 2 * np.pi
        elif lon - lon1 > np.pi:
            lon -= 2 * np.pi

        points.append((np.degrees(lat), np.degrees(lon)))

    return points

# 羽田空港の緯度経度
haneda_lat = 35.5533  # haneda_lat = 33.4484（Phoenixの緯度）
haneda_lon = 139.7811 # haneda_lon = -112.0740（Phoenixの経度）
# Dallasの緯度経度
dallas_lat = 32.7767
dallas_lon = -96.7970

# 地図を作成
m = folium.Map(location=[haneda_lat, haneda_lon], zoom_start=3)

# 羽田空港にマーカーを追加
folium.Marker([haneda_lat, haneda_lon], popup='Haneda Airport, Tokyo').add_to(m)

# 大円上の点を計算
num_points = 100
great_circle = great_circle_points((haneda_lat, haneda_lon), (dallas_lat, dallas_lon), num_points)

# 大円を地図に追加
folium.PolyLine(locations=great_circle, color='blue').add_to(m)

# Dallasにマーカーを追加（修正後の位置に変更）
folium.Marker([great_circle[-1][0], great_circle[-1][1]], popup='Dallas, TX').add_to(m)

# 地図を表示
m

# 地図をHTMLファイルとして保存
# m.save("flight_map.html")