WTOAAM/dataAnalysisBusiness/demo/SCADA_10min_category_1.py

import os
import re
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator#设定固定刻度

def scada_10min_category():
    turbine_number=24
    
    fpath = 'D:/赵雅丽/实习/算法/min_scada_LuoTuoGou/72/'
    # 定义一个正则表达式来匹配纯数字文件名且扩展名为.csv  
    pattern = re.compile(r'^\d+\.csv$')
    
    # 列出指定路径下的所有文件和文件夹  
    files_in_dir = os.listdir(fpath)
    for file in files_in_dir:  
        # 使用正则表达式匹配文件名  
        if pattern.match(file):  
            # 拼接文件的完整路径  
            fname = os.path.join(fpath, file)
            # 读取csv文件，保持原始变量名而不忽略任何行
            scada_10min = pd.read_csv(fname)
        
            # 显示数据        
            time_stamp = scada_10min.loc[:,['时间']] #dataframe
            active_power = scada_10min.loc[:,['变频器电网侧有功功率']]
            wind_speed = scada_10min.loc[:,['风速']]
            LM = pd.concat([time_stamp,active_power,wind_speed],axis=1)  #dataframe
            # lm=LM.values #array
    
            xx = data_label(LM,fpath)#dataframe
            mergedTable = pd.concat([scada_10min,xx],axis=1)#合并dataframe
            D = mergedTable[mergedTable['lab'] == 1]#选择为1的行
            ws = D["风速"]#series
            ap = D["变频器电网侧有功功率"]
            ##绘图
            # fig = plt.figure(figsize=(10,6),dpi=500)  #figsize是图形大小，dpi像素
            plt.scatter(ws,ap,s=8,c='black',marker='.',label='好点')
            # x_major_locator=MultipleLocator(5)
            # y_major_locator=MultipleLocator(500)
            # ax=plt.gca()
            # ax.xaxis.set_major_locator(x_major_locator)
            # ax.yaxis.set_major_locator(y_major_locator)
            # plt.xlim((0,30))
            # plt.ylim((0,2200))
            # plt.tick_params(labelsize=20)
            # # plt.grid(c='dimgray',alpha=0.2)
            # plt.xlabel("V/(m$·$s$^{-1}$)",fontsize=20)
            # plt.ylabel("P/kW",fontsize=20)

            # # plt.savefig(r'D:\赵雅丽\研究生学习资料\学习资料\劣化度健康度\spyder\大论文\图\风速-功率.jpg',bbox_inches='tight')
            # plt.show()
        
        
        
def data_label(x1,x2):   # LM:T P V  path:文件获取路径
    fpath2 = x2
    fname2 = os.path.join(fpath2, "info.csv") #读取数据文件2（额定风速额定功率等）
    # 参数na_filter=False仅阻止了pandas自动检测这些缺失值，并不能忽略  
    # 但请注意，pandas没有直接的'omitrow'选项，如果需要忽略包含缺失值的行，需要在后续处理中处理
    turbine_info = pd.read_csv(fname2, na_filter=False)
    # 删除包含任何缺失值的行  
    turbine_info = turbine_info.dropna() 
    
    PRated = turbine_info.loc[:,["额定功率"]] #dataframe
    VCutOut = turbine_info.loc[:,["切出风速"]]  
    VCutIn = turbine_info.loc[:,["切入风速"]]  
    VRated = turbine_info.loc[:,["额定风速"]]
    
    #网格法确定风速风向分区数量，功率方向分区数量
    Labeled_March809 = x1
    APower = Labeled_March809["active_power"]  #series读入有功功率
    WSpeed = Labeled_March809["wind_speed"]  #读入风速
    maxP=np.max(APower)
    intervalP=25  #ceil(PRated*0.01)#功率分区间隔为额定功率的1%
    intervalwindspeed=0.25  #风速分区间隔0.25m/s
    #初始化
    PNum = 0  
    TopP = 0   
    # 根据条件计算PNum和TopP  
    if maxP >= PRated:  
        PNum = math.floor(maxP / intervalP) + 1  
        TopP = math.floor((maxP - PRated) / intervalP) + 1  
    else:  
        PNum = math.floor(PRated / intervalP)  
        TopP = 0   
    VNum = math.ceil(VCutOut / intervalwindspeed)  
    SM1 = Labeled_March809.shape  
    AA1 = SM1[0]  #运行数据的条数
    lab = [[0] for _ in range(AA1)]  #创建全0空列表
    lab = pd.DataFrame(lab,columns=['lab'])
    Labeled_March809 = pd.concat([Labeled_March809,lab],axis=1)  #在tpv后加一列标签列
    SM = Labeled_March809.shape #(52561,4)
    AA = SM[0]  
    #存储功率大于0的运行数据
    #标识功率为0的点，标识-1
    DzMarch809_0 = np.zeros(AA, 3)  #array(52561,3)
    nCounter1 = 1
    Point_line=np.zeros(AA,1)
    #考虑到很多功率小于10的数据存在，将<10的功率视为0
    for i in range(AA):
        if (APower[i] > 10) & (WSpeed[i] > 0):
            nCounter1 += 1   #共有nCounter1个功率大于0的正常数据
            DzMarch809_0[nCounter1-1, 0] = WSpeed[i]  
            DzMarch809_0[nCounter1-1, 1] = APower[i]  
            Point_line[nCounter1-1] = i+1  # 记录nCounter1记下的数据在原始数据中的位置  
        if APower[i] <= 10: 
            Labeled_March809[i,SM[1]-1] = -1  # 功率为0标识为-1  array类型
    # 截取DzMarch809_0中实际存储的数据  其他全为0
    DzMarch809 = DzMarch809_0[:nCounter1, :]  
    #统计各网格落入的散点个数
    XBoxNumber = np.ones((PNum, VNum),dtype=int)  #(86 100)
    nWhichP = 0
    nWhichV = 0

    # 循环遍历DzMarch809中的有效数据  
    for i in range(nCounter1):  
        
        # 查找功率所在的区间  
        for m in range(1, PNum + 1):  # 注意Python的range是左闭右开的，所以需要+1  
            if (DzMarch809[i,1] > (m - 1) * intervalP) and (DzMarch809[i,1] <= m * intervalP):  
                nWhichP = m  
                break  
          
        # 查找风速所在的区间  
        for n in range(1, VNum + 1):  # 同样需要+1  
            if (DzMarch809[i, 0] > (n - 1)*intervalwindspeed) and (DzMarch809[i, 0] <= n*intervalwindspeed):  
                nWhichV = n  
                break  
          
        # 如果功率和风速都在有效区间内，增加对应网格的计数  
        if (nWhichP > 0) and (nWhichV > 0):  
            XBoxNumber[nWhichP - 1, nWhichV - 1] += 1  # 注意Python的索引是从0开始的，所以需要减1  
    # XBoxNumber现在包含了每个网格的计数[PNum行, VNum列]

    for m in range(1,PNum+1):
        for n in range(1,VNum+1):
            XBoxNumber[m-1,n-1] = XBoxNumber[m-1,n-1] - 1

    #在功率方向将网格内散点绝对个数转换为相对百分比，备用
    PBoxPercent = np.zeros((PNum, VNum))  #(86 100) #计算后会出现浮点型，所以不能定义int类型
    PBinSum = np.zeros((PNum,1),dtype=int)
    for i in range(1,PNum+1):
        for m in range(1,VNum+1):
            PBinSum[i-1] = PBinSum[i-1] + XBoxNumber[i-1,m-1] 
        for m in range(1,VNum+1):
            if PBinSum[i-1]>0:
                PBoxPercent[i-1,m-1] = (XBoxNumber[i-1,m-1] / PBinSum[i-1])*100
    #在风速方向将网格内散点绝对个数转换为相对百分比，备用          
    VBoxPercent = np.zeros((PNum, VNum))  #(86 100) #计算后会出现浮点型，所以不能定义int类型
    VBinSum = np.zeros((VNum,1),dtype=int)
    for i in range(1,VNum+1):
        for m in range(1,PNum+1):
            VBinSum[i-1] = VBinSum[i-1] + XBoxNumber[m-1,i-1] 
        for m in range(1,PNum+1):
            if VBinSum[i-1]>0:
                VBoxPercent[m-1,i-1] = (XBoxNumber[m-1,i-1] / VBinSum[i-1])*100
    # VBoxPercent PBoxPercent 左上-右下
    # 将数据颠倒一下  左下-右上         第一行换为倒数第一行 方便可视化
    InvXBoxNumber = np.zeros((PNum,VNum),dtype = int)
    InvPBoxPercent = np.zeros((PNum,VNum),dtype = float)
    InvVBoxPercent = np.zeros((PNum,VNum),dtype = float)
    for m in range(1,PNum+1):
        for n in range(1,VNum+1):
            InvXBoxNumber[m-1,n-1] = XBoxNumber[PNum-(m-1)-1,n-1]
            InvPBoxPercent[m-1,n-1] = PBoxPercent[PNum-(m-1)-1,n-1]
            InvVBoxPercent[m-1,n-1] = VBoxPercent[PNum-(m-1)-1,n-1]
    
    #以水平功率带方向为准，分析每个水平功率带中，功率主带中心，即找百分比最大的网格位置。
    PBoxMaxIndex = np.zeros((PNum,1),dtype = int)  #水平功率带最大网格位置索引
    PBoxMaxP = np.zeros((PNum,1),dtype = float)       #水平功率带最大网格百分比
    for m in range(1,PNum+1):
        PBoxMaxIndex[m-1] = np.argmax(PBoxPercent[m-1, :])   #argmax返回最大值的索引
        PBoxMaxP[m-1] = np.max(PBoxPercent[m-1, :])
    #以垂直风速方向为准，分析每个垂直风速带中，功率主带中心，即找百分比最大的网格位置。
    VBoxMaxIndex = np.zeros((VNum,1),dtype = int)  
    VBoxMaxV = np.zeros((VNum,1),dtype = float)       
    for m in range(1,VNum+1):
        VBoxMaxIndex[m-1] = np.argmax(VBoxPercent[:, m-1])   
        VBoxMaxV[m-1] = np.max(VBoxPercent[:, m-1])
    
    #切入风速特殊处理，如果切入风速过于偏右，向左拉回
    if PBoxMaxIndex[0]>14:                     #第一个值对应的是风速最小处 即切入风速
        PBoxMaxIndex[0] = 9 
    #以水平功率带方向为基准，进行分析
    DotDense = np.zeros(PNum)   #每一水平功率带的功率主带包含的网格数
    DotDenseLeftRight = np.zeros((PNum,2))  #存储每一水平功率带的功率主带以最大网格为中心，向左，向右扩展的网格数
    DotValve = 90  #从中心向左右对称扩展网格的散点百分比和的阈值。
    PDotDenseSum = 0
    for i in range(PNum - TopP):  # 从最下层水平功率带开始，向上分析到特定的功率带  
        PDotDenseSum = PBoxMaxP[i]  # 以中心最大水平功率带为基准，向左向右对称扩展网格，累加各网格散点百分比  
        iSpreadRight = 1  
        iSpreadLeft = 1  
          
        while PDotDenseSum < DotValve:  
            if (PBoxMaxIndex[i] + iSpreadRight) < VNum-1-1:  
                PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] + iSpreadRight]  # 向右侧扩展  
                iSpreadRight += 1  
            else:
                break  
              
            if (PBoxMaxIndex[i] - iSpreadLeft) > 0:  
                PDotDenseSum += PBoxPercent[i, PBoxMaxIndex[i] - iSpreadLeft]  # 向左侧扩展  
                iSpreadLeft += 1  
            else:  
                break  
        iSpreadRight = iSpreadRight-1
        iSpreadLeft = iSpreadLeft-1
        #向左右扩展完毕
        DotDenseLeftRight[i, 0] = iSpreadLeft  # 左  
        DotDenseLeftRight[i, 1] = iSpreadRight  # 右  
        DotDense[i] = iSpreadLeft + iSpreadRight + 1  # 记录向左向右扩展的个数及每个功率仓内网格的个数  
    # 此时DotDense和DotDenseLeftRight数组已经包含了所需信息    
    #各行功率主带右侧宽度的中位数最具有代表性（因为先右后左）
    DotDenseWidthLeft = np.zeros((PNum-TopP))
    for i in range(PNum-TopP):
        DotDenseWidthLeft[i] = DotDenseLeftRight[i,1]  #DotDenseLeftRight[i,1]:向右延伸个数
    MainBandRight = np.median(DotDenseWidthLeft) #计算中位数
    
    # 初始化变量  
    PowerLimit = np.zeros(PNum, dtype=int)  # 各水平功率带是否为限功率标识，1：是；0：不是  
    WidthAverage = 0  # 功率主带右侧平均宽度  
    WidthAverage_L = 0  # 功率主带左侧平均宽度  
    WidthVar = 0  # 功率主带方差（此变量在提供的代码中并未使用）  
    PowerLimitValve = 6  # 限功率主带判别阈值  
    N_Pcount = 20  # 阈值  
      
    nCounterLimit = 0  # 限功率的个数  
    nCounter = 0  # 正常水平功率带的个数  
      
    # 循环遍历水平功率带，从第1个到第PNum-TopP个  
    for i in range(PNum - TopP):  
        # 如果向右扩展网格数大于阈值，且该水平功率带点总数大于20，则标记为限功率带  
        if (DotDenseLeftRight[i, 1] > PowerLimitValve) and (PBinSum[i] > N_Pcount):  
            PowerLimit[i] = 1  
            nCounterLimit += 1  #限功率的个数
          
        # 如果向右扩展网格数小于等于阈值，则累加右侧宽度（左侧宽度在代码中似乎有误）  
        if DotDenseLeftRight[i, 1] <= PowerLimitValve:  
            WidthAverage += DotDenseLeftRight[i, 1]  # 统计正常水平功率带右侧宽度
            WidthAverage_L += DotDenseLeftRight[i,1]   #统计正常水平功率带左侧宽度
            nCounter += 1  
    # 计算平均宽度  
    WidthAverage /= nCounter if nCounter > 0 else 1  # 避免除以0的情况  
    WidthAverage_L /= nCounter if nCounter > 0 else 1   

    #计算正常（即非限功率）水平功率带的功率主带宽度的方差，以此来反映从下到上宽度是否一致
    WidthVar = 0  # 功率主带宽度的方差   
    for i in range(PNum - TopP):  
        # 如果向右扩展网格数小于等于阈值，则计算当前宽度与平均宽度的差值平方  
        if DotDenseLeftRight[i, 1] <= PowerLimitValve:  
            WidthVar += (DotDenseLeftRight[i, 1] - WidthAverage) ** 2  
    # 计算方差（注意：除以nCounter-1是为了得到样本方差）  
    WidthVar = np.sqrt(WidthVar / (nCounter - 1) if nCounter > 1 else 0)  # 避免除以0的情况

    #各水平功率带，功率主带的风速范围，右侧扩展网格数*2*0.25
    PowerBandWidth = WidthAverage*intervalwindspeed+WidthAverage_L*intervalwindspeed

    # 对限负荷水平功率带的最大网格进行修正  
    for i in range(1, PNum - TopP+1):  
        if (PowerLimit[i] == 1) and (abs(PBoxMaxIndex[i] - PBoxMaxIndex[i - 1]) > 5):  
            PBoxMaxIndex[i] = PBoxMaxIndex[i - 1] + 1  
      
    # 输出各层功率主带的左右边界网格索引  
    DotDenseInverse = np.flipud(DotDenseLeftRight)  # 上下翻转数组以得到反向顺序  
      
    # 计算功率主带的左右边界  
    CurveWidthR = np.ceil(WidthAverage) + 2  # 功率主带的右边界 + 2  
    CurveWidthL = np.ceil(WidthAverage_L) + 2  # 功率主带的左边界 + 2  
      
    # 网格是否为限功率网格的标识数组  
    BBoxLimit = np.zeros((PNum, VNum), dtype=int)  
    # 标记限功率网格  
    for i in range(2, PNum - TopP):  
        if PowerLimit[i] == 1:
            BBoxLimit[i, int(PBoxMaxIndex[i] + CurveWidthR + 1):VNum] = 1

    # 初始化数据异常需要剔除的网格标识数组  
    BBoxRemove = np.zeros((PNum, VNum), dtype=int)  
    # 标记需要剔除的网格  
    for m in range(PNum - TopP): 
        for n in range(int(PBoxMaxIndex[m]) + int(CurveWidthR), VNum):  # 注意Python中的索引从0开始，因此需要减去1  
            BBoxRemove[m, n] = 1 
        # 功率主带左侧的超发网格，从最大索引向左直到第一个网格  
        
        for n in range(int(PBoxMaxIndex[m]) - int(CurveWidthL)+1, 0, -1):  # 使用range的步长参数来实现从右向左的迭代  
            BBoxRemove[m, n-1] = 2  # 注意Python中的索引从0开始，因此需要减去1

    # 初始化变量  
    CurveTop = np.zeros((2, 1), dtype=int)  
    CurveTopValve = 1  # 网格的百分比阈值  
    BTopFind = 0  
    mm = 0  
    #确定功率主带的左上拐点，即额定风速位置的网格索引
    CurveTop = np.zeros((2, 1), dtype=int)  
    CurveTopValve = 1  # 网格的百分比阈值  
    BTopFind = 0  
    mm = 0   
    for m in range(PNum - TopP, 0, -1):  # 注意Python的range是左闭右开区间，所以这里从PNum-TopP开始到1（不包括0）  
        for n in range(int(np.floor(int(VCutIn) / intervalwindspeed)), VNum - 1):  # 使用floor函数来向下取整  
            if (VBoxPercent[m, n - 1] < VBoxPercent[m, n]) and (VBoxPercent[m, n] <= VBoxPercent[m, n + 1]) and (XBoxNumber[m, n] >= 3):   
                CurveTop[0] = m  
                CurveTop[1] = n  #[第80个,第40个]
                BTopFind = 1  
                mm = m  # mm是拐点所在功率仓，对应其index
                break  # 找到后退出内层循环  
        if BTopFind == 1:  
            break  # 找到后退出外层循环
            
    IsolateValve = 3  #功率主带右侧孤立点占比功率仓阈值 3%
    # 遍历功率仓和网格  
    for m in range(PNum - TopP):    
        for n in range(int(PBoxMaxIndex[m]) + int(CurveWidthR), VNum):  
            # 检查PBoxPercent是否小于阈值，如果是，则标记BBoxRemove为1  
            if PBoxPercent[m, n] < IsolateValve:   
                BBoxRemove[m, n] = 1
    #功率主带顶部宽度
    CurveWidthT = np.floor((maxP - PRated) / intervalP) + 1  
    # 标记额定功率以上的超发点(PNum-PTop之间)  
    for m in range(PNum - TopP, PNum):   
        for n in range(VNum):  
            BBoxRemove[m, n] = 3
      
    # 标记功率主带拐点左侧的欠发网格  
    for m in range(mm-1, PNum - TopP): 
        for n in range(int(CurveTop[1]) - 2):
            BBoxRemove[m, n] = 2    # BBoxRemove数组现在包含了根据条件标记的超发点和欠发网格的信息

    #以网格的标识，决定该网格内数据的标识。
    # DzMarch809Sel数组现在包含了每个数据点的标识
    DzMarch809Sel = np.zeros(nCounter1, dtype=int)  # 初始化标识数组   
    nWhichP = 0  
    nWhichV = 0  
    nBadA = 0   
    for i in range(nCounter1):  
        for m in range( PNum ):   
            if (DzMarch809[i, 1] > m * intervalP) and (DzMarch809[i, 1] <= (m+1) * intervalP):  
                nWhichP = m  #m记录的是index
                break  
        for n in range( VNum ):  # 注意Python的range是左闭右开区间，所以这里到VNum+1  
            if DzMarch809[i, 0] > ((n+1) * intervalwindspeed - intervalwindspeed/2) and DzMarch809[i, 0] <= ((n+1) * intervalwindspeed + intervalwindspeed / 2):  
                nWhichV = n  #index
                break  
        if nWhichP >= 0 and nWhichV >= 0:  
            if BBoxRemove[nWhichP, nWhichV] == 1:   
                DzMarch809Sel[i] = 1  
                nBadA += 1  
            elif BBoxRemove[nWhichP, nWhichV] == 2:  
                DzMarch809Sel[i] = 2  
            elif BBoxRemove[nWhichP , nWhichV] == 3:  
                DzMarch809Sel[i] = 0  # 额定风速以上的超发功率点认为是正常点，不再标识  
    # DzMarch809Sel数组现在包含了每个数据点的标识

    ##############################滑动窗口方法
    # 存储限负荷数据  
    PVLimit = np.zeros((nCounter1, 3))  #存储限负荷数据  %第3列用于存储限电的点所在的行数
    nLimitTotal = 0  
    nWindowLength = 6   #滑动窗口长度设置为6
    LimitWindow = np.zeros(nWindowLength)  #滑动窗口空列表
    UpLimit = 0    #上限
    LowLimit = 0   #下限
    PowerStd = 30  # 功率波动方差  
    nWindowNum = np.floor(nCounter1/nWindowLength) #6587
    PowerLimitUp = PRated - 100  
    PowerLimitLow = 100  

    # 循环遍历每个窗口  
    for i in range(int(nWindowNum)):  
        start_idx = i * nWindowLength  
        end_idx = start_idx + nWindowLength  
        LimitWindow = DzMarch809[start_idx:end_idx, 1]  # 提取当前窗口的数据  
          
        # 检查窗口内所有数据是否在功率范围内  
        bAllInAreas = np.all(LimitWindow >= PowerLimitLow) and np.all(LimitWindow <= PowerLimitUp)  
        if not bAllInAreas:  
            continue  
          
        # 计算方差上下限  
        UpLimit = LimitWindow[0] + PowerStd  
        LowLimit = LimitWindow[0] - PowerStd  
          
        # 检查窗口内数据是否在方差范围内  
        bAllInUpLow = np.all(LimitWindow >= LowLimit) and np.all(LimitWindow <= UpLimit)  
        if bAllInUpLow:  
            # 标识窗口内的数据为限负荷数据  
            DzMarch809Sel[start_idx:end_idx] = 4  
              
            # 存储限负荷数据  
            for j in range(nWindowLength):  
                PVLimit[nLimitTotal, :2] = DzMarch809[start_idx + j, :2]  
                PVLimit[nLimitTotal, 2] = Point_line[start_idx + j]  # 对数据进行标识  
                nLimitTotal += 1  
    # PVLimit现在包含了限负荷数据，nLimitTotal是限负荷数据的总数
    
    
    #将功率滑动窗口主带平滑化
    # 初始化锯齿平滑的计数器  
    nSmooth = 0  
    # 遍历除了最后 TopP+1 个元素之外的所有 PBoxMaxIndex 元素  
    for i in range(PNum - TopP - 1):  
        PVLeftDown = np.zeros(2)  
        PVRightUp = np.zeros(2)  
        # 检查当前与下一个 PBoxMaxIndex 之间的距离是否大于等于1  
        if PBoxMaxIndex[i + 1] - PBoxMaxIndex[i] >= 1:  
            # 计算左下和右上顶点的坐标  
            PVLeftDown[0] = (PBoxMaxIndex[i]+1 + CurveWidthR) * 0.25 - 0.125  
            PVLeftDown[1] = (i) * 25  
            PVRightUp[0] = (PBoxMaxIndex[i+1]+1 + CurveWidthR) * 0.25 - 0.125  
            PVRightUp[1] = (i+1) * 25  
              
            # 遍历 DzMarch809 数组  
            for m in range(nCounter1):  
                # 检查当前点是否在锯齿区域内  
                if (DzMarch809[m, 0] > PVLeftDown[0]) and (DzMarch809[m, 0] < PVRightUp[0]) and (DzMarch809[m, 1] > PVLeftDown[1]) and (DzMarch809[m, 1] < PVRightUp[1]):
                    # 检查斜率是否大于对角连线  
                    if (DzMarch809[m, 1] - PVLeftDown[1]) / (DzMarch809[m, 0] - PVLeftDown[0]) > (PVRightUp[1] - PVLeftDown[1]) / (PVRightUp[0] - PVLeftDown[0]):
                        # 如果在锯齿左上三角形中，则选中并增加锯齿平滑计数器  
                        DzMarch809Sel[m] = 0  
                        nSmooth += 1  
    # DzMarch809Sel 数组现在＋包含了锯齿平滑的选择结果，nSmooth 是选中的点数
    ###################################存储数据
    # 存储好点  
    nCounterPV = 0  # 初始化计数器  
    PVDot = np.zeros((nCounter1, 3))  # 初始化存储好点的数组  nCounter1是p>0的数
    for i in range(nCounter1):  
        if DzMarch809Sel[i] == 0:  
            nCounterPV += 1  
            PVDot[nCounterPV-1, :2] = DzMarch809[i, :2]  
            PVDot[nCounterPV-1, 2] = Point_line[i]  # 好点 Point_line记录nCounter1在原始数据中的位置 
    nCounterVP = nCounterPV  
     
    # 对所有数据中的好点进行标注    
    for i in range(nCounterVP):  
        Labeled_March809[int(PVDot[i, 2] - 1), (SM[1]-1)] = 1  # 注意Python的索引从0开始，并且需要转换为整数索引  
     
    # 存储坏点  
    nCounterBad = 0  # 初始化计数器  
    PVBad = np.zeros((nCounter1, 3))  # 初始化存储坏点的数组  
    for i in range(nCounter1):  
        if DzMarch809Sel[i] in [1, 2, 3]:  
            nCounterBad += 1  
            PVBad[nCounterBad-1, :2] = DzMarch809[i, :2]  
            PVBad[nCounterBad-1, 2] = Point_line[i]  
        
    # 对所有数据中的坏点进行标注  
    for i in range(nCounterBad):  
        Labeled_March809[int(PVBad[i, 2] - 1),(SM[1]-1)] = 5  # 坏点标识  

    # 对所有数据中的限电点进行标注   
    for i in range(nLimitTotal):  
        Labeled_March809[int(PVLimit[i, 2] - 1),(SM[1]-1)] = 4  # 限电点标识  

    # 对所有的数据点进行标注  
    # Labeled_March809是array，提取所第四列的值保存为dataframe
    A = Labeled_March809[:,3]
    A=pd.DataFrame(A,columns=['lab'])
    return A


# scada_10min_category()