如何用Python绘制跟随飞行路径的3D动态带状图？

1. 数据读取与解析 #

日志每行采用键值对格式，用正则表达式提取字段最直接：

import re
import numpy as np

def parse_log_file(file_path):
    data = {
        'timestamp': [],
        'pitch': [],
        'roll': [],
        'altitude_ft': [],
        'latitude': [],
        'longitude': [],
        'heading': [],
        'turn_rate': []
    }
    # 匹配日志行的正则模式
    pattern = re.compile(r'Time: (\d+\.\d+) alt \(ft\): (\d+\.\d+) lat\(deg\): (\-?\d+\.\d+) lon\(deg\): (\-?\d+\.\d+) pitch\(deg\): (\-?\d+\.\d+) roll\(deg\): (\-?\d+\.\d+) heading\(deg\): (\-?\d+\.\d+) turnspeed: (\-?\d+\.\d+)')
    
    with open(file_path, 'r') as f:
        for line in f:
            match = pattern.match(line.strip())
            if match:
                data['timestamp'].append(float(match.group(1)))
                data['altitude_ft'].append(float(match.group(2)))
                data['latitude'].append(float(match.group(3)))
                data['longitude'].append(float(match.group(4)))
                data['pitch'].append(float(match.group(5)))
                data['roll'].append(float(match.group(6)))
                data['heading'].append(float(match.group(7)))
                data['turn_rate'].append(float(match.group(8)))
    
    # 转numpy数组方便后续计算
    for key in data:
        data[key] = np.array(data[key])
    return data

2. 地理坐标转3D惯性坐标系 #

将经纬度、海拔转换为笛卡尔坐标（简化为球体模型，若需高精度可改用WGS84椭球公式）：

def latlonalt_to_cartesian(lat_deg, lon_deg, alt_ft):
    # 地球平均半径（米）
    EARTH_RADIUS = 6371000
    # 英尺转米
    alt_m = alt_ft * 0.3048
    lat_rad = np.radians(lat_deg)
    lon_rad = np.radians(lon_deg)
    
    x = (EARTH_RADIUS + alt_m) * np.cos(lat_rad) * np.cos(lon_rad)
    y = (EARTH_RADIUS + alt_m) * np.cos(lat_rad) * np.sin(lon_rad)
    z = (EARTH_RADIUS + alt_m) * np.sin(lat_rad)
    return x, y, z

调用转换函数处理数据：

flight_data = parse_log_file('flight_log.txt')
x, y, z = latlonalt_to_cartesian(flight_data['latitude'], flight_data['longitude'], flight_data['altitude_ft'])

3. 生成带状路径几何数据 #

根据飞行方向（heading）生成路径两侧的带状端点，假设带状宽度为100米：

def get_flight_direction(heading_deg):
    heading_rad = np.radians(heading_deg)
    # 正北对应Y轴正方向，东对应X轴正方向，计算归一化飞行方向向量
    dx = np.sin(heading_rad)
    dy = np.cos(heading_rad)
    dz = 0  # 基础版忽略俯仰影响，如需加入可结合pitch调整
    norm = np.sqrt(dx**2 + dy**2 + dz**2)
    return dx/norm, dy/norm, dz/norm

band_width = 100  # 带状宽度（米）
left_x, left_y, left_z = [], [], []
right_x, right_y, right_z = [], [], []

for i in range(len(x)):
    dx, dy, dz = get_flight_direction(flight_data['heading'][i])
    # 计算垂直于飞行方向的向量（地面垂直方向叉乘飞行方向）
    up_vec = np.array([0, 0, 1])
    cross_vec = np.cross(np.array([dx, dy, dz]), up_vec)
    cross_vec = cross_vec / np.linalg.norm(cross_vec) * (band_width / 2)
    
    # 记录左右端点坐标
    left_x.append(x[i] + cross_vec[0])
    left_y.append(y[i] + cross_vec[1])
    left_z.append(z[i] + cross_vec[2])
    right_x.append(x[i] - cross_vec[0])
    right_y.append(y[i] - cross_vec[1])
    right_z.append(z[i] - cross_vec[2])

# 转numpy数组
left_x, left_y, left_z = np.array(left_x), np.array(left_y), np.array(left_z)
right_x, right_y, right_z = np.array(right_x), np.array(right_y), np.array(right_z)

若需体现roll姿态的影响，可通过旋转矩阵将垂直向量按roll角旋转，让带状随飞机倾斜。

4. 动态3D绘图 #

使用matplotlib的animation模块实现逐帧更新的动态效果：

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 初始化3D画布
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel('X (m)')
ax.set_ylabel('Y (m)')
ax.set_zlabel('Z (m)')
ax.set_title('3D Dynamic Flight Path Band')

# 设置坐标轴范围（根据实际数据调整）
ax.set_xlim(np.min(x)-1000, np.max(x)+1000)
ax.set_ylim(np.min(y)-1000, np.max(y)+1000)
ax.set_zlim(np.min(z)-1000, np.max(z)+1000)

# 初始化绘图对象
path_line, = ax.plot([], [], [], 'b-', linewidth=2)
band_left, = ax.plot([], [], [], 'r-', linewidth=1)
band_right, = ax.plot([], [], [], 'r-', linewidth=1)
fill_poly = ax.plot_surface([], [], [], color='blue', alpha=0.3)

def update(frame):
    # 更新中心路径线
    path_line.set_data(x[:frame+1], y[:frame+1])
    path_line.set_3d_properties(z[:frame+1])
    
    # 更新带状左右边界
    band_left.set_data(left_x[:frame+1], left_y[:frame+1])
    band_left.set_3d_properties(left_z[:frame+1])
    band_right.set_data(right_x[:frame+1], right_y[:frame+1])
    band_right.set_3d_properties(right_z[:frame+1])
    
    # 更新填充的带状面
    global fill_poly
    fill_poly.remove()
    X = np.array([left_x[:frame+1], right_x[:frame+1]])
    Y = np.array([left_y[:frame+1], right_y[:frame+1]])
    Z = np.array([left_z[:frame+1], right_z[:frame+1]])
    fill_poly = ax.plot_surface(X, Y, Z, color='blue', alpha=0.3)
    
    return path_line, band_left, band_right, fill_poly

# 创建动画，interval为帧间隔（毫秒）
ani = FuncAnimation(fig, update, frames=len(x), interval=50, blit=True)

plt.show()

5. 优化方向 #

• 动态宽度调整：可根据turn rate实时修改带状宽度，增强机动动作的视觉表现
• 姿态融合：引入旋转矩阵，结合roll/pitch参数调整带状的倾斜角度
• 性能优化：对大数据量做降采样处理，减少帧数量提升流畅度
• 坐标系适配：若惯性坐标系定义不同，可修改坐标转换公式调整轴方向

内容的提问来源于stack exchange，提问作者spacestar