import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist

# Step 1: 从 CSV 文件读取数据
# 假设文件为 'flight_data.csv'，文件中包含 flight_id, x, y, z, t, heading, v
data = pd.read_csv('../data.csv')  # 读取 CSV 文件

# Step 2: 按 flight_id 分组，并且对每个轨迹按时间 t 排序
# 获取唯一的 flight_id
flight_ids = data['flight_id'].unique()

# 存储不同轨迹的列表
trajectories = {}

# 对每个 flight_id 进行数据分组，并按时间 t 排序
for flight_id in flight_ids:
    flight_data = data[data['flight_id'] == flight_id]  # 获取该 flight_id 的所有数据
    flight_data_sorted = flight_data.sort_values(by='t')  # 按时间 t 排序
    trajectories[flight_id] = flight_data_sorted

# Step 2: 获取轨迹特征（保持每条轨迹的多维坐标）
trajectory_features = []

for flight_id, trajectory in trajectories.items():
    trajectory_points = trajectory[['x', 'y']].values  # 获取轨迹的所有 (x, y) 点
    trajectory_features.append(trajectory_points)  # 保持多维坐标，不展平

# 将所有轨迹点转换为一个 NumPy 数组（每条轨迹是一个二维数组）
trajectory_features = np.array(trajectory_features)

# Step 3: 计算轨迹之间的欧式距离
# 这里使用 cdist 计算轨迹之间的距离，我们需要将每条轨迹转换为相同长度的形状
# 设定轨迹点数量相同，这里假设每条轨迹有 N 个点

# 假设轨迹的最大点数是 N（例如，最大100个点）
N =200  # 这里可以根据实际数据调整
distances = np.zeros((len(trajectories), len(trajectories)))

# 对每对轨迹，计算轨迹之间的欧式距离
for i in range(len(trajectory_features)):
    for j in range(i+1, len(trajectory_features)):
        trajectory_i = trajectory_features[i]
        trajectory_j = trajectory_features[j]
        
        # 如果轨迹长度不同，可以进行补齐或截断，这里采用截断策略
        # if len(trajectory_i) > N:
        #     trajectory_i = trajectory_i[:N]
        # if len(trajectory_j) > N:
        #     trajectory_j = trajectory_j[:N]
        
        # 计算轨迹之间的欧式距离
        dist = np.sum(np.sqrt(np.sum((trajectory_i - trajectory_j) ** 2, axis=1)))  # 按逐点距离求和
        distances[i, j] = dist
        distances[j, i] = dist
        # distances.append(dist)
from sklearn.metrics import silhouette_score
import warnings
warnings.filterwarnings("ignore")
# Step 3: 手肘法计算
inertia = []  # 用于存储不同 K 值的 inertia（误差平方和）
K_range = range(2, 11)  # 聚类数 K 从 1 到 10
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(distances)
    inertia.append(kmeans.inertia_)

# 绘制手肘法图
plt.figure(figsize=(8, 6))
plt.plot(K_range, inertia, marker='o')
plt.title('Elbow Method for Optimal K')
plt.xlabel('Number of Clusters (K)')
plt.ylabel('Inertia')
plt.show()

# Step 4: 计算轮廓系数
best_k = 4  # 假设从手肘法中选择的 K 值是 3
kmeans = KMeans(n_clusters=best_k, random_state=42)
kmeans_labels = kmeans.fit_predict(distances)

# 计算轮廓系数
sil_score = silhouette_score(distances, kmeans_labels)
print(f"Silhouette Score for K={best_k}: {sil_score:.4f}")

# Step 5: 展示聚类结果
plt.figure(figsize=(8, 6))
colors = plt.cm.get_cmap('tab10', best_k)

for i, flight_id in enumerate(trajectories):
    trajectory = trajectories[flight_id]
    plt.plot(trajectory['x'], trajectory['y'], color=colors(kmeans_labels[i]), label=f'Cluster {kmeans_labels[i]+1}' if i == 0 else "")

plt.title('KMeans Clustering of Flight Trajectories')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.legend()
plt.show()

# 输出聚类结果
for flight_id, label in zip(trajectories.keys(), kmeans_labels):
    print(f"Flight {flight_id}: Cluster {label + 1}")
我想加入z标签这一列的数据，并将绘图结果改为xyz三维形式，可以帮我输出修改后的代码吗