import numpy as np
import pandas as pd
from sklearn.cluster import KMeans, OPTICS
import matplotlib.pyplot as plt
from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score
from scipy.spatial.distance import directed_hausdorff
import warnings

# 读取数据
data = pd.read_csv(r'C:\Users\29244\Desktop\data.csv')
flight_ids = data['flight_id'].unique()
trajectories = {}

for flight_id in flight_ids:
    flight_data = data[data['flight_id'] == flight_id]
    flight_data_sorted = flight_data.sort_values(by='t')
    trajectories[flight_id] = flight_data_sorted

trajectory_features = []

for flight_id, trajectory in trajectories.items():
    trajectory_points = trajectory[['x', 'y']].values
    trajectory_features.append(trajectory_points)

trajectory_features = np.array(trajectory_features)

# 初始化距离矩阵
N = 200
distances_euclidean = np.zeros((len(trajectories), len(trajectories)))
distances_directed_hausdorff = np.zeros((len(trajectories), len(trajectories)))

# 欧式距离计算
def euclidean_distance(series1, series2):
    return np.linalg.norm(series1 - series2)

# 单向豪斯多夫距离计算
def directed_hausdorff_distance(series1, series2):
    return directed_hausdorff(series1, series2)[0]  

# 计算距离矩阵
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_euclidean = euclidean_distance(trajectory_i, trajectory_j)
        dist_directed_hausdorff = directed_hausdorff_distance(trajectory_i, trajectory_j)

        distances_euclidean[i, j] = dist_euclidean
        distances_euclidean[j, i] = dist_euclidean
        distances_directed_hausdorff[i, j] = dist_directed_hausdorff
        distances_directed_hausdorff[j, i] = dist_directed_hausdorff

# KMeans聚类
def perform_kmeans_clustering(distances, k):
    model = KMeans(n_clusters=k, random_state=42)
    labels = model.fit_predict(distances)
    return labels

# OPTICS聚类
def perform_optics_clustering(distances):
    model = OPTICS(min_samples=2, metric='precomputed')
    labels = model.fit_predict(distances)
    return labels

k_values = [4, 4]
results = {}

# 聚类并绘制结果图
for dist, name in zip([distances_euclidean, distances_directed_hausdorff], ['欧式距离', '单向豪斯多夫距离']):
    for method in ['KMeans', 'OPTICS']:
        if method == 'KMeans':
            labels = perform_kmeans_clustering(dist, k_values[0])
        else:
            labels = perform_optics_clustering(dist)

        silhouette = silhouette_score(dist, labels) if method == 'KMeans' else '未定义'
        ch_index = calinski_harabasz_score(dist, labels) if method == 'KMeans' else '未定义'
        db_index = davies_bouldin_score(dist, labels) if method == 'KMeans' else '未定义'
        
        results[(name, method)] = {
            '轮廓系数': silhouette,
            'CH指数': ch_index,
            'DB指数': db_index
        }
        
        colors = plt.cm.get_cmap('tab10', k_values[0]) if method == 'KMeans' else plt.cm.get_cmap('tab10', len(set(labels)) - (1 if -1 in labels else 0))
        plt.figure(figsize=(8, 6))
        
        for i, flight_id in enumerate(trajectories):
            trajectory = trajectories[flight_id]
            plt.plot(trajectory['x'], trajectory['y'], color=colors(labels[i]), label=f'聚类 {labels[i] + 1}' if i == 0 else "")
        
        plt.title(f'{name}聚类飞行轨迹（{method}）', fontproperties='SimHei')
        plt.xlabel('经度', fontproperties='SimHei')
        plt.ylabel('纬度', fontproperties='SimHei')
        plt.xlim([data['x'].min(), data['x'].max()])
        plt.ylim([data['y'].min(), data['y'].max()])
        plt.legend(prop={'family': 'SimHei'})  
        plt.show()

# 输出结果指标
for key, metrics in results.items():
    print(f"{key[0]} - {key[1]}:")
    print(f"  轮廓系数: {metrics['轮廓系数']:.4f}" if metrics['轮廓系数'] != '未定义' else "  轮廓系数: 未定义")
    print(f"  CH指数: {metrics['CH指数']:.4f}" if metrics['CH指数'] != '未定义' else "  CH指数: 未定义")
    print(f"  DB指数: {metrics['DB指数']:.4f}" if metrics['DB指数'] != '未定义' else "  DB指数: 未定义\n")
以上是一个机场一天中进场航班的二维航迹聚类，下面我给你评价系数结果，研究结果形成实验报告，分析轨迹聚类结果，讨论不同相似性度量以及聚类方法的优缺点欧式距离 - KMeans:
  轮廓系数: 0.6487
  CH指数: 805.0958
  DB指数: 0.4674
欧式距离 - OPTICS:
  轮廓系数: 未定义
  CH指数: 未定义
  DB指数: 未定义

单向豪斯多夫距离 - KMeans:
  轮廓系数: 0.4993
  CH指数: 291.3701
  DB指数: 0.7714
单向豪斯多夫距离 - OPTICS:
  轮廓系数: 未定义
  CH指数: 未定义
  DB指数: 未定义