# !usr/bin/env python # -*- coding:utf-8 _*- """ @Author:Lijiaxing @File:data_analysis.py @Time:2023/4/24 15:16 """ import os.path import pandas as pd # from mpl_toolkits.basemap import Basemap from scipy.signal import savgol_filter import numpy as np import matplotlib.pyplot as plt from scipy.cluster.hierarchy import dendrogram, linkage, fcluster from sklearn.metrics import silhouette_samples, silhouette_score def paint_others(y): """ 绘制其他数据 """ plt.plot([j for j in range(y)], y) # 添加标题和标签 plt.xlabel('x') plt.ylabel('y') # 显示图形 plt.show() def compute_cos_similarity(a, b): """ 计算两个向量的余弦相似度 :param a: 向量a :param b: 向量b :return: 余弦相似度值 """ return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b)) def compute_pearsonr(a): """ 计算数据皮尔逊相关系数并返回相似度矩阵 :param a: 数据格式为n*m的矩阵，n为数据个数，m为数据维度 :return: 返回相似度矩阵,数据格式为n*n的矩阵 """ return np.corrcoef(a) def compute_distance(a, b): """ 计算两个向量的欧式距离 :param a: :param b: :return: 返回两个向量的欧式距离 """ return np.linalg.norm(a - b) def hierarchical_clustering(data, threshold, similarity_func): """ 层次聚类,使用工具包scipy.cluster.hierarchy中的linkage和fcluster函数进行层次聚类 :param data: 二维数据，格式为n*m的矩阵，n为数据个数，m为数据维度 :param threshold: 阈值，当两个数据的距离小于阈值时，将两个数据归为一类，阈值为根据相似度矩阵层次聚类后的类别距离阈值，可根据需求进行调整，可大于1 :param similarity_func: 相似度计算函数，用于计算两个数据的相似度，可以进行替换，若替换为计算距离的函数需对内部进行修改 :return: 返回聚类结果，格式为n*1的矩阵，n为数据个数，每个数据的值为该数据所属的类别 """ # 计算数据的相似度矩阵 similarity_matrix = similarity_func(data) # 计算数据的距离矩阵 distance_matrix = 1 - similarity_matrix # 进行层次聚类返回聚类结果 Z = linkage(distance_matrix, method='ward') # 根据相似度阈值获取聚类结果 clusters = fcluster(Z, t=threshold, criterion='distance') # 画出层次聚类树形结构 fig = plt.figure(figsize=(5, 3)) dn = dendrogram(Z) plt.show() # clusters[42] = 1 silhouette = silhouette_samples(np.abs(distance_matrix), clusters, metric='euclidean') silhouette1 = silhouette_score(np.abs(distance_matrix), clusters, metric='euclidean') print(f"平均轮廓系数为：{silhouette1}, 单个样本的轮廓系数：{silhouette}") return clusters class DataAnalysis: """ 数据分析类 """ def __init__(self, data_length, data_start, data_end): """ 初始化 :param data_length: 分析数据段长度 :param data_start: 分析数据段开始位置 :param data_end: 分析数据段结束位置 """ # 原始风机功率数据傅里叶变换滤波后的数据 self.ori_turbine_fft = None # 原始风机功率数据片段 self.ori_turbine_pic = None # 聚类结果 self.cluster = None # 风机功率差分平滑后的结果 self.smooth_turbine_diff = None # 风机功率差分变化情况 self.diff_change = None # 风机功率差分 self.turbine_diff = None # 全部风机数据 self.turbine = None # 风机的标号顺序 self.turbine_id = list(range(102, 162)) # b1b4 = [142, 143, 144, 145] # self.turbine_id = [id for id in self.turbine_id if id not in b1b4] # 风机功率数据15分钟级别 self.power_15min = None # 风机经纬度信息 self.info = None # 使用数据长度 self.data_length = data_length # 使用数据开始位置 self.data_start = data_start # 使用数据结束位置 self.data_end = data_end # 导入数据 self.load_data(normalize=True) # 计算风机功率差分 self.compute_turbine_diff() def load_data(self, normalize=False): """ 加载数据 :return: """ self.info = pd.read_csv('../data-process/data/风机信息.csv', encoding='utf-8') # power_15min = pd.read_csv('../data/power_15min.csv') # for i in range(len(power_15min)): # if power_15min.loc[i, 'C_REAL_VALUE'] == -9999: # # 方便在曲线中看出缺失数据位置 # power_15min.loc[i, 'C_REAL_VALUE'] = -34.56789 # self.power_15min = power_15min turbine_path = '../data-process/data/output_filtered_csv_files/turbine-{}.csv' self.turbine, turbines = {}, [] for i in self.turbine_id: self.turbine[i] = pd.read_csv(turbine_path.format(i))[20:].reset_index(drop=True) if normalize is True: self.normalize() def normalize(self): turbines = [self.turbine[i].values[:, 1:].astype(np.float32) for i in self.turbine_id] turbines = np.vstack(turbines) mean, std = np.mean(turbines, axis=0), np.std(turbines, axis=0) for i in self.turbine_id: c_time = self.turbine[i]['C_TIME'] self.turbine[i] = (self.turbine[i].iloc[:, 1:] - mean) / std self.turbine[i].insert(loc=0, column='C_TIME', value=c_time) return self.turbine def compute_turbine_diff(self): """ 计算风机功率差分 :return: """ turbine_diff = [] ori_turbine_pic = [] for turbine_i in self.turbine_id: ori = np.array(self.turbine[turbine_i]['C_WS'].values[self.data_start:self.data_end + 1]) diff_array = np.diff(ori) smoothness_value = np.std(diff_array) print("turbine-{}的平滑度是：{}".format(turbine_i, round(smoothness_value, 2))) turbine_diff.append(diff_array) ori_turbine_pic.append(self.turbine[turbine_i]['C_WS'].values[self.data_start:self.data_end]) self.ori_turbine_pic = ori_turbine_pic self.turbine_diff = turbine_diff diff_change = [] for diff_i in turbine_diff: single_diff_change = [] for diff_i_i in diff_i: if diff_i_i > 0: single_diff_change.append(1) elif diff_i_i < 0: single_diff_change.append(-1) else: single_diff_change.append(0) diff_change.append(single_diff_change) self.diff_change = diff_change self.ori_turbine_fft = [self.turbine_fft(i + 1) for i in range(len(self.ori_turbine_pic))] # 平滑 self.turbine_smooth(window_size=21) def paint_map(self): """ 绘制经纬度地图 :return: """ lats = self.info['纬度'].values lons = self.info['经度'].values map = Basemap() # 绘制海岸线和国家边界 map.drawcoastlines() map.drawcountries() # 绘制经纬度坐标 map.drawmeridians(range(0, 360, 30)) map.drawparallels(range(-90, 90, 30)) # 绘制点 x, y = map(lons, lats) map.plot(x, y, 'bo', markersize=10) # 显示图表 plt.show() def paint_power15min(self): """ 绘制15分钟功率曲线 :return: """ plt.plot(self.power_15min['C_REAL_VALUE']) # 设置图表标题和轴标签 plt.title('Data Time Change Curve') plt.xlabel('Date') plt.ylabel('Value') # 显示图表 plt.show() def paint_lats_lons(self): """ 绘制经纬度图 :return: """ x = self.info['纬度'].values y = self.info['经度'].values # 绘制散点图 fig, ax = plt.subplots() plt.scatter(x, y) for i, txt in enumerate(self.info['id'].values): ax.annotate(txt, (x[i], y[i])) # 设置图表标题和轴标签 plt.xlabel('lats') plt.ylabel('lons') # 显示图表 plt.show() def similarity_score(self, turbine_diff, threshold=0.5): """ 使用余弦相似度计算相似度分数并返回相似度大于阈值的index矩阵 :param turbine_diff: 需要计算相似的矩阵，数据格式n*m,n为数据条数，m为数据维数 :param threshold: 相似度阈值 :return: 返回相似计算后的矩阵 """ similarity = {i: [] for i in range(49)} similarity_index = {i: [] for i in range(49)} for turbine_i in range(49): for turbine_j in range(49): cos_similarity = compute_cos_similarity(turbine_diff[turbine_i], turbine_diff[turbine_j]) similarity[turbine_i].append(cos_similarity) if cos_similarity > threshold: similarity_index[turbine_i].append(turbine_j) return similarity_index def paint_turbine(self, paint_default=True): """ 绘制风机地理位置图 :param paint_default:默认True，绘制聚类后每个类别的数据折线图 :return: None """ # y = self.info['纬度'].values # x = self.info['经度'].values # # fig, ax = plt.subplots(figsize=(15, 15)) # # plt.scatter(x, y, c=self.cluster) # for i, txt in enumerate(self.info['C_ID'].values): # ax.annotate(txt, (x[i], y[i])) # 设置图表标题和轴标签 # plt.xlabel('lons') # plt.ylabel('lats') # plt.legend() # # # 显示图表 # plt.savefig('analysis_img/turbine_cluster.png') # plt.show() plt.figure(figsize=(20, 10)) cmap = plt.get_cmap('viridis') linestyle= ['solid', 'dashed'] for i in range(max(self.cluster)): cluster, cluster_fft = [], [] for j, item in enumerate(self.cluster): if item == i + 1: cluster.append(self.ori_turbine_pic[j]) cluster_fft.append(self.ori_turbine_fft[j]) cluster_fft = np.average(cluster_fft, axis=0) cluster = np.average(cluster, axis=0) diff_array = np.diff(cluster) smoothness_value = np.std(diff_array) print("聚类-{}的平滑度是：{}".format(i+1, smoothness_value)) color = cmap(i*200) plt.subplot(2, 1, 1) plt.plot([j for j in range(len(cluster))], cluster, color=color, label='cluster'+str(i), linestyle=linestyle[i]) plt.subplot(2, 1, 2) plt.plot([j for j in range(len(cluster_fft))], cluster_fft, color=color, label='cluster'+str(i), linestyle=linestyle[i]) # 添加图例 plt.legend() # 显示图形 plt.savefig('analysis_img/cluster/clusters.png') plt.show() if paint_default: for i in range(max(self.cluster)): self.paint_turbine_k(i + 1) # 画出聚类中每个风机的曲线 def turbine_smooth(self, window_size=50): """ 使用滑动平均平滑数据。 参数： data -- 需要平滑的数据，numpy数组类型 window_size -- 滑动窗口大小，整数类型 返回值： smooth_data -- 平滑后的数据，numpy数组类型 """ # weights = np.repeat(1.0, window_size) / window_size smooth_data = [] for turbine_diff_i in self.turbine_diff: smooth_y = savgol_filter(turbine_diff_i, window_length=window_size, polyorder=3) smooth_data.append(smooth_y) # smooth_data.append(np.convolve(turbine_diff_i, weights, 'valid')) self.smooth_turbine_diff = smooth_data def paint_turbine_k(self, k): """ 绘制第k聚类的风机数据折线图 :param k: :return: """ pic_label = [] y = [] plt.figure(figsize=(20, 10)) cmap = plt.get_cmap('viridis') for i, item in enumerate(self.cluster): if item == k: pic_label.append('turbine-' + str(self.turbine_id[i])) y.append(self.ori_turbine_fft[i]) for i in range(len(y)): color = cmap(i / 10) plt.plot([j for j in range(len(y[i]))], y[i], color=color, label=pic_label[i]) # 添加标签和标题 plt.xlabel('x') plt.ylabel('y') plt.title('Cluster {}'.format(k)) # 添加图例 plt.legend() # 显示图形 plt.savefig('analysis_img/cluster/cluster_{}.png'.format(k)) plt.show() def turbine_fft(self, k): """ 对第k台原始风机数据进行傅里叶变换，并绘制变换前后曲线 :param k: 数据读入时的风机顺序index，从1开始 :return: 傅里叶变换清洗后的数据，数据格式 """ y = self.ori_turbine_pic t = np.linspace(0, 1, self.data_length) signal = y[k - 1] # 进行傅里叶变换 freq = np.fft.fftfreq(len(signal), t[1] - t[0]) spectrum = np.fft.fft(signal) spectrum_abs = np.abs(spectrum) threshold = np.percentile(spectrum_abs, 98) indices = spectrum_abs > threshold spectrum_clean = indices * spectrum # 进行傅里叶逆变换 signal_clean = np.fft.ifft(spectrum_clean) # plt.figure(figsize=(20, 10)) # # # 绘制时域信号 # plt.subplot(4, 1, 1) # plt.plot(t, signal) # plt.title(self.turbine_id[k-1]) # # # 绘制频域信号 # plt.subplot(4, 1, 2) # plt.plot(freq, np.abs(spectrum)) # # # 绘制过滤后的频域信号 # plt.subplot(4, 1, 3) # plt.plot(freq, np.abs(spectrum_clean)) # # # 绘制过滤后的时域信号 # plt.subplot(4, 1, 4) # plt.plot(t, signal_clean) # # plt.savefig('analysis_img/fft/{}_turbine_fft.png'.format(self.turbine_id[k-1])) # plt.show() return signal_clean def paint_double(self, i, j): """ 绘制两台风机的数据变换对比 :param i: 风机数据读入时数据编号，从1开始 :param j: 风机数据读入时数据编号，从1开始 :return: """ y = self.ori_turbine_fft x = [index for index in range(self.data_length)] data_i = y[i - 1] data_j = y[j - 1] plt.figure(figsize=(20, 10)) plt.plot(x, data_i, label='turbine {}'.format(self.turbine_id[i - 1]), linestyle='solid') plt.plot(x, data_j, label='turbine {}'.format(self.turbine_id[j - 1]), linestyle='dashed') plt.title('{} and {}'.format(i, j)) plt.legend() plt.savefig('analysis_img/{}_{}_turbine.png'.format(self.turbine_id[i - 1], self.turbine_id[j - 1])) plt.show() def process_ori_data(self): """ 对原始数据进行处理，聚类和绘图 :return: """ self.turbine_clusters(self.ori_turbine_fft) self.paint_turbine() def turbine_clusters(self, data=None): """ 风机数据聚类，聚类信息保存在self.cluster中 :param data: 默认为空，也可以使用其他数据聚类，并体现在绘图中， 数据格式:二维数据n*m，n为数据条数，m为每条数据维数 :return: None """ if data is None: cluster = hierarchical_clustering(self.turbine_diff, threshold=1.4, similarity_func=compute_pearsonr) # 层次聚类 else: cluster = hierarchical_clustering(data, threshold=0.8, similarity_func=compute_pearsonr) self.cluster = cluster # 在这里保存cluster变量 from cluster_analysis import cluster_power_list_file, cluster_power_list_folder output_path = '../data-process/data/cluster_power/' cluster_power_list_file(self.cluster, self.turbine_id, input_path='../data-process/data/output_filtered_csv_files/', output_path=output_path) cluster_power_list_folder(self.cluster, self.turbine_id, input_path='../data-process/data/continuous_data/', output_path=output_path) data_analysis = DataAnalysis(data_length=9773, data_start=0, data_end=9773) data_analysis.process_ori_data() data_analysis.paint_double(1, 56)