WTOAAM/dataAnalysisBusiness/algorithm/yawErrorDensityAnalyst.py

import os
import pandas as pd
import numpy as np
import plotly.graph_objects as go
from algorithmContract.confBusiness import *
from algorithmContract.contract import Contract
from behavior.analystWithGoodBadLimitPoint import AnalystWithGoodBadLimitPoint
from scipy.stats import binned_statistic_2d
from scipy.stats import skew, kurtosis
from utils.jsonUtil import JsonUtil
from scipy.stats import norm, gaussian_kde

class YawErrorDensityAnalyst(AnalystWithGoodBadLimitPoint):
    """
    风电机组动态偏航策略分析
    """

    def typeAnalyst(self):
        return "yaw_error_density"

    def selectColumns(self):
        return [Field_DeviceCode,Field_Time,Field_WindSpeed,Field_ActiverPower,Field_YawError]

    def turbinesAnalysis(self,  outputAnalysisDir, conf: Contract, turbineCodes):
        dictionary = self.processTurbineData(turbineCodes,conf,self.selectColumns())
        dataFrameMerge = self.userDataFrame(dictionary,conf.dataContract.configAnalysis,self)
        turbineInfos = self.common.getTurbineInfos(conf.dataContract.dataFilter.powerFarmID, turbineCodes, self.turbineInfo)
        results=self.yawErrorAnalysis(dataFrameMerge, turbineInfos,outputAnalysisDir, conf)
        return results

    def yawErrorAnalysis(self, dataFrameMerge: pd.DataFrame, turbineModelInfo: pd.Series, outputAnalysisDir, conf: Contract):
        # 检查所需列是否存在
        required_columns = {Field_ActiverPower, Field_YawError,Field_WindSpeed}
        if not required_columns.issubset(dataFrameMerge.columns):
            raise ValueError(f"DataFrame缺少必要的列。需要的列有: {required_columns}")

        result_rows = []
        grouped = dataFrameMerge.groupby(
            [Field_NameOfTurbine, Field_CodeOfTurbine])
                # 定义固定的颜色映射列表
        fixed_colors = [
                "#3E409C",
                "#476CB9",
                "#3586BF",
                "#4FA4B5",
                "#52A3AE",
                "#60C5A3",
                "#85D0AE",
                "#A8DCA2",
                "#CFEE9E",
                "#E4F39E",
                "#EEF9A7",
                "#FBFFBE",
                "#FDF1A9",
                "#FFE286",
                "#FFC475",
                "#FCB06C",
                "#F78F4F",
                "#F96F4A",
                "#E4574C",
                "#CA3756",
                "#AF254F"
        ]

        # 将 fixed_colors 转换为 Plotly 的 colorscale 格式
        fixed_colorscale = [
            [i / (len(fixed_colors) - 1), color] for i, color in enumerate(fixed_colors)
        ]

        for name, group in grouped:
            dataFrame=group[Field_YawError].abs() <= 45
            # yawerror=np.mean(dataFrame[Field_YawError])

            df = self.calculateYawError(group)
            
            # 这里的 dropna 会修改 df
            df.dropna(inplace=True)
            
            # 1. 之前加的防护：如果是空的，跳过
            if df.empty:
                # print(f"Warning: Turbine {name[0]} has no valid data after screening, skip the analysis.")
                continue
            
            counts = df['density'].value_counts()
            count_0 = counts.get(0, 0)  # 获取 density 为 0 的数量，如果没有 0 则返回 0
            count_1 = counts.get(1, 0)  # 获取 density 为 1 的数量，如果没有 1 则返回 0

            # 2. 这里的筛选可能再次导致 df 变空
            df = df[df["density"]>0]

            # 【重要】再次检查是否为空，或者数据量太少
            # KDE 至少需要 2 个点才能计算，且方差不能为 0
            if df.empty or len(df) < 2:
                 # print(f"Warning: Turbine {name[0]} has insufficient data for density plot.")
                 continue

            mean_X = np.mean(df["x"])
            std_X = np.std(df["x"])
            mediann_X= np.median(df["x"])
            skewness_X = skew(df["x"])
            kurtosis_X = kurtosis(df["x"])
            max_X = np.max(df["x"])
            min_X = np.min(df["x"])
            result={
                'mean_X':[mean_X],
                'std_X': [std_X],
                'mediann_X':[mediann_X],
                'skewness_X':[skewness_X],
                'kurtosis_X':[kurtosis_X],
                'max_X':[max_X],
                'min_X':[min_X]
            }
            result = pd.DataFrame(result)

            # 用密度作为颜色绘制散点图
            fig = go.Figure()
            fig.add_trace(go.Scattergl(
                x=df["x"],
                y=df["y"],
                mode='markers',
                marker=dict(
                    size=3,
                    opacity=0.7,
                    color=df["density"],
                    colorscale=fixed_colorscale,
                    showscale=True,
                )
            ))

            fig.update_layout(
                xaxis_title='对风角度',
                yaxis_title='风速',
                title=f'动态偏航误差分析-{name[0]}',
                xaxis=dict(range=[-20, 20]),  # 限制横坐标范围为 -20 到 20
                yaxis=dict(range=[0, 25])
            )

            # 确保从 Series 中提取的是具体的值
            engineTypeCode = turbineModelInfo.get(Field_MillTypeCode, "")
            if isinstance(engineTypeCode, pd.Series):
                engineTypeCode = engineTypeCode.iloc[0]

            engineTypeName = turbineModelInfo.get(Field_MachineTypeCode, "")
            if isinstance(engineTypeName, pd.Series):
                engineTypeName = engineTypeName.iloc[0]
            # 构建最终的JSON对象
            json_output = {
                "analysisTypeCode": "动态偏航误差",
                "engineCode":  engineTypeCode,
                "engineTypeName": engineTypeName,
                "xaixs": "对风角度(°)",
                "yaixs": "风速(m/s)",
                "data": [{
                "engineName": name[0],
                "engineCode": name[1],
                "title":f'动态偏航误差分析-{name[0]}',
                "xData": df["x"].tolist(),
                "yData": df["y"].tolist(),
                "colorbar": df["density"].tolist(),
                "colorbartitle": "密度"

            }]
            }

            x_range = np.linspace(-30, 30, 1000)  # 先定义 x 轴
            try:
                # 尝试计算 KDE
                # 要求：数据点 > 1 且 方差 != 0
                if df["x"].nunique() > 1: 
                    kde = gaussian_kde(df["x"])
                    pdf_data = kde(x_range)
                else:
                    # 如果所有数据都一样（方差为0），KDE 无法计算
                    # 直接给 0 或者抛出异常进入 except
                    raise ValueError("Data variance is zero")
            except Exception as e:
                # print(f"KDE calculation failed for {name[0]}: {e}. Using zero distribution.")
                # 如果失败，生成全 0 的平直线，保证 JSON 结构完整
                pdf_data = np.zeros_like(x_range)

            # 构建最终的JSON对象2
            json_output2 = {
                "analysisTypeCode": "动态偏航误差",
                "engineCode":  engineTypeCode,
                "engineTypeName": engineTypeName,
                "xaixs": "对风角度(°)",
                "yaixs": "概率密度函数",
                "data": [{
                "engineName": name[0],
                "engineCode": name[1],
                "title":f'概率密度函数-{name[0]}',
                "xData": x_range.tolist(),
                "yData": pdf_data.tolist(),
                "xrange":[-30,30]

            }]
            }

            # Save to file
            # filePathOfImage = os.path.join(outputAnalysisDir, f"{name[0]}.png")
            # fig.write_image(filePathOfImage, scale=3)
            # filePathOfHtml = os.path.join(outputAnalysisDir, f"{name[0]}.html")
            # fig.write_html(filePathOfHtml)


            # 将JSON对象保存到文件np.mean(dataFrame[Field_YawError])
            output_json_path = os.path.join(outputAnalysisDir, f"{name[0]}.json")
            with open(output_json_path, 'w', encoding='utf-8') as f:
                import json
                json.dump(json_output, f, ensure_ascii=False, indent=4)

            # 将JSON对象2保存到文件
            output_json_path2 = os.path.join(outputAnalysisDir, f"PDF-{name[0]}.json")
            with open(output_json_path2, 'w', encoding='utf-8') as f:
                import json
                json.dump(json_output2, f, ensure_ascii=False, indent=4)

            # 如果需要返回DataFrame，可以包含文件路径
            result_rows.append({
                Field_Return_TypeAnalyst: self.typeAnalyst(),
                Field_PowerFarmCode: conf.dataContract.dataFilter.powerFarmID,
                Field_Return_BatchCode: conf.dataContract.dataFilter.dataBatchNum,
                Field_CodeOfTurbine: name[1],
                Field_Return_FilePath: output_json_path,
                Field_Return_IsSaveDatabase: True
            })
            result_rows.append({
                Field_Return_TypeAnalyst: self.typeAnalyst(),
                Field_PowerFarmCode: conf.dataContract.dataFilter.powerFarmID,
                Field_Return_BatchCode: conf.dataContract.dataFilter.dataBatchNum,
                Field_CodeOfTurbine: name[1],
                Field_Return_FilePath: output_json_path2,
                Field_Return_IsSaveDatabase: True
            })

            # result_rows.append({
            #     Field_Return_TypeAnalyst: self.typeAnalyst(),
            #     Field_PowerFarmCode: conf.dataContract.dataFilter.powerFarmID,
            #     Field_Return_BatchCode: conf.dataContract.dataFilter.dataBatchNum,
            #     Field_CodeOfTurbine: name[1],
            #     Field_Return_FilePath: filePathOfImage,
            #     Field_Return_IsSaveDatabase: False
            # })

            # result_rows.append({
            #     Field_Return_TypeAnalyst: self.typeAnalyst(),
            #     Field_PowerFarmCode: conf.dataContract.dataFilter.powerFarmID,
            #     Field_Return_BatchCode: conf.dataContract.dataFilter.dataBatchNum,
            #     Field_CodeOfTurbine: name[1],
            #     Field_Return_FilePath: filePathOfHtml,
            #     Field_Return_IsSaveDatabase: True


            # })

        result_df = pd.DataFrame(result_rows)

        return result_df

    def calculateYawError(self, dataFrame: pd.DataFrame):

        dataFrame = dataFrame.dropna(
            subset=[Field_NameOfTurbine, Field_YawError, Field_ActiverPower,Field_WindSpeed])

        filtered_dataFrame = dataFrame[(dataFrame[Field_YawError].abs() <= 30)&(dataFrame[Field_WindSpeed] >= 0)&(dataFrame[Field_WindSpeed] <= 25)]
        
        # 如果筛选结果为空，直接返回 None，不再进行后续计算，防止 binned_statistic_2d 报错
        if filtered_dataFrame.empty:
            return pd.DataFrame()
        
        x=filtered_dataFrame[Field_YawError]
        y=filtered_dataFrame[Field_WindSpeed]
        # data = np.column_stack((x, y))  # 合并为两列数组
        # 使用 binned_statistic_2d 来计算散点的密度分布
        binSize_x = 60
        binSize_y = 50
        counts, x_edges, y_edges, binnumber = binned_statistic_2d(x, y, values=None, statistic='count', bins=[binSize_x, binSize_y])

        # 将数据密度转化为颜色值
        binX = np.digitize(x, x_edges) - 1
        binY = np.digitize(y, y_edges) - 1

        # 删除超出范围的下标
        validIdx = (binX >= 0) & (binX < binSize_x) & (binY >= 0) & (binY < binSize_y)

        # 获取有效下标的密度
        density = np.zeros(len(x))
        density[validIdx] = counts[binX[validIdx], binY[validIdx]]
        # 将结果保存到 result 中
        result = {
            'x': x,
            'y': y,
            'density': density
        }

        # 将 result 转换为 DataFrame
        result_df = pd.DataFrame(result)

        return result_df