温室气体数据订证之PICARRO（python、bz2file、pandas、matplotlib)

大家好，我是讯享网，很高兴认识大家。

今年其实还是写了不少代码的，但来去还是在玩pandas，就今年重点解决的一个问题来分享。

数据说明

温室气体设备PICARRO打包的数据格式为bz2文件，包括了CO2,CH4,CO，水汽，切换的阀口等数据。

数据处理的总设计

数据的读入部分

使用bz2file包来读bz2压缩文件，返回dataframe格式的数据，部分代码如下。

#Author Wu Dongqiao 2020.11.18 import os import numpy as np import pandas as pd import time import matplotlib.pyplot as plt from matplotlib.font_manager import FontProperties import bz2file font = FontProperties(fname=r"C:\\WINDOWS\\Fonts\\simsun.ttc", size=14) import warnings warnings.filterwarnings("ignore") #1.读入数据模块 def RDT(path): #读co2和co的数据函数 def getdata_CO2CO(path): file = bz2file.open(path, "r") line = file.readline() TIME, FLAG, CO2, CO = [], [], [], [] while True: line = file.readline() if not line: # 读到文件尾后退出 break # 科学字符串转科学直数用eval再转十进制用%d time = str(line.split()[0], 'utf-8') + ' ' + str(line.split()[1], 'utf-8') flag = int('%d' % eval(str(line.split()[11], 'utf-8'))) co2 = float('%f' % eval(str(line.split()[15], 'utf-8'))) co = float('%f' % eval(str(line.split()[16], 'utf-8'))) TIME.append(time) FLAG.append(flag) CO2.append(co2) CO.append(co) df = pd.DataFrame({ 
   'time': TIME, 'flag': FLAG, 'CO2': CO2, 'CO': CO}) df['time'] = pd.to_datetime(df['time'], format='%m/%d/%y %H:%M:%S.%f') df = df.set_index('time') return df #读CH4数据的函数 def getdata_CH4(path): file = bz2file.open(path, "r") line = file.readline() TIME, CH4 = [], [] while True: line = file.readline() if not line: break # 科学字符串转科学直数用eval再转十进制用%d time = str(line.split()[0], 'utf-8') + ' ' + str(line.split()[1], 'utf-8') # flag=int('%d'%eval(str(line.split()[11], 'utf-8'))) ch4 = float('%f' % eval(str(line.split()[16], 'utf-8'))) TIME.append(time) CH4.append(ch4) df = pd.DataFrame({ 
   'time': TIME, 'CH4': CH4}) df['time'] = pd.to_datetime(df['time'], format='%m/%d/%y %H:%M:%S.%f') df = df.set_index('time') df.dropna() return df # 设置了空的dataframe用于保存数据 data_CO2CO = pd.DataFrame(columns=('time', 'flag', 'CO2', 'CO')) data_CO2CO = data_CO2CO.set_index('time') data_CH4 = pd.DataFrame(columns=('time', 'CH4')) data_CH4 = data_CH4.set_index('time') for root, dirs, files in os.walk(path): for tempFile in files: print(tempFile) path = os.path.join(root, tempFile) if path[-16] == 'C': data_CO2CO = pd.concat([data_CO2CO, getdata_CO2CO(path)], axis=0) elif path[-16] == 'H': data_CH4 = pd.concat([data_CH4, getdata_CH4(path)], axis=0) # 合并了CO2,CO,CH4的数据在一个df data = pd.concat([data_CO2CO, data_CH4], axis=1) return data

讯享网

时间处理模块部分，时间去重，处理切换阀口后的两分钟左右的数据

在这里插入图片描述

每次换阀后带来一分钟左右的混气段，需删除。

讯享网#2.时间处理模块,时间去重，处理切换后的120秒左右的数据，缺失处理 def del_c(data): temptime=data.index.tolist() for ii in range(len(data)-1): if (temptime[ii + 1] - temptime[ii]).seconds > 200: data.loc[temptime[ii+1]:temptime[ii+120]]=np.nan data=data.dropna() return data def TMK(data): #分类成样气段和标气断，并分别去掉转换的前120个混合数据） data1=data[(data['flag'] == 1) | (data['flag']==2)] data4=data[(data['flag'] == 4) | (data['flag']==5)] data_12=del_c(data1) data_45 =del_c(data4) return data_12,data_45

异常值处理模块

异常值处理，根据要素制定的合理性范围来去除异常值，实际上还应该根据水气、压力、流量等状态来进行处理，这里只放范围处理部分

def SDT(df): # 处理因换冷凝管或其它原因导致的数据异常,范围参考 df.loc[(df['CO2'] > 480) | (df['CO2'] < 380)] = np.nan df.loc[(df['CO'] > 2) | (df['CO'] < 0.1)] = np.nan df.loc[(df['CH4'] > 5) | (df['CH4'] < 1.8)] = np.nan df.dropna() return df

线性校准模块

订证最重要的部分，使用y=ax+b外标法定正，使用两种已知浓度标气来进行标定，首先已知浓度的标气浓度作为一个常量。

讯享网 '''1.W和T标气值，两点校准的参考，应该从TXT里读取，每次更换标气需填加或修改''' CRDS_W_CO2 = 402.19 CRDS_W_CH4 = 1.9414 CRDS_W_CO = 0.3901 CRDS_T_CO2 = 398.39 CRDS_T_CH4 = 1.898 CRDS_T_CO = 0.1821

校准模块

def CONC(data_SAM,data_STA): ''''2.利用原数据，读取标气段，每段一个平均值（去掉前60秒）''' S=data_STA temptime=S.index.tolist() Time,S_CO2, S_CO, S_CH4,S_flag= [],[],[],[],[] #新表的列 x, y, z ,n= 0, 0, 0,0 #过程变量 for ii in range(len(S)-1): #.index.tolist n=n+1 x=x+S.iloc[ii]['CO2'] y=y+S.iloc[ii]['CO'] z=z+S.iloc[ii]['CH4'] if (temptime[ii+1]-temptime[ii]).seconds>1200:#当下个时次的标定时间大于上次的20min作为分段标准 Time.append(temptime[ii]) # print(temptime[ii]) S_CO2.append(np.round(x/n,2)) S_CO.append(np.round(y / n, 4)) S_CH4.append(np.round(z / n, 4)) S_flag.append(S.iloc[ii-2]['flag']) #部分端口为5口，最后几个flag显示为4.再往上倒数2个 x, y, z, n = 0, 0, 0, 0 S_data=pd.DataFrame({ 
   'time':Time,'CO2':S_CO2,'CO':S_CO,'CH4':S_CH4,'flag':S_flag}) S_data=S_data.set_index('time') S_data.to_csv('S_data.csv') #保存新生成的校准表 '''3.按顺序历遍校准表，每个值传出开始时间（上一段结束时间），结束时间（本段时间），计算校准参数a, b 历遍校准表，找到原数据表对应时间段，按y = ax + b处理数据''' #样气数据不分层 data_sample=data_SAM # data_sample.to_csv('data.csv') # 求系数斜率公式 def asc(x1, y1, x2, y2): a = round((y1 - y2) / (x1 - x2), 4) return a # 求截矩公式 def inte(x1, y1, x2, y2): b = round((x1 * y2 - x2 * y1) / (x1 - x2), 4) return b # 订式公式 def lin(a, b, x): y = a * x + b return y #历遍并根据时段分别进行校准 temptime=S_data.index.tolist() for ii in range(len(S_data)-1): time_S=temptime[ii] #校准的开始时间 time_E=temptime[ii+1] #校准的结束时间 print(time_S) if S_data.iloc[ii]['flag']==4: data_WS = S_data.iloc[ii] data_TS = S_data.iloc[ii+1] if S_data.iloc[ii]['flag']==5: data_TS=S_data.iloc[ii] data_WS = S_data.iloc[ii+1] # print(S_data.iloc[ii]['flag']) # 求系数（截矩） CO2_b = inte(CRDS_W_CO2, data_WS['CO2'], CRDS_T_CO2, data_TS['CO2']) CO_b = inte(CRDS_W_CO, data_WS['CO'], CRDS_T_CO, data_TS['CO']) CH4_b = inte(CRDS_W_CH4, data_WS['CH4'], CRDS_T_CH4, data_TS['CH4']) # 求系数斜率 CO2_a = asc(CRDS_W_CO2, data_WS['CO2'], CRDS_T_CO2, data_TS['CO2']) CO_a = asc(CRDS_W_CO, data_WS['CO'], CRDS_T_CO, data_TS['CO']) CH4_a = asc(CRDS_W_CH4, data_WS['CH4'], CRDS_T_CH4, data_TS['CH4']) #订正 data_sample.loc[time_S:time_E]['CO2']=data_sample.loc[time_S:time_E]['CO2'].apply(lambda x:lin(CO2_a, CO2_b, x)) data_sample.loc[time_S:time_E]['CO']=data_sample.loc[time_S:time_E]['CO'].apply(lambda x:lin(CO_a, CO_b, x)) data_sample.loc[time_S:time_E]['CH4']=data_sample.loc[time_S:time_E]['CH4'].apply(lambda x:lin(CH4_a, CH4_b, x)) # data_sample.to_csv('DZ_data.csv') # print(data_sample) return data_sample

主程序调用

讯享网#Author Wu Dongqiao 2020.12.2 if __name__=="__main__": #指定数据文件夹，每个文件夹最好只放少于一个月的文件，秒数据量比较大，每次处理所使用的时间比较多。 path=r'D:\data\picarro\' '''1.调用读入数据模块''' data=RDT(path) '''2.时间处理，去掉转换后的120秒数据''' data_SAM,data_STA=TMK(data) '''3.调用异常值处理模块''' data_SAM=SDT(data_SAM) data_STA=SDT(data_STA) #分为高低层数据后，5min平均 data_sampleH = data_SAM[data_SAM['flag'] == 1].resample('5min').mean() #标气一分钟平均 data_STA_4=data_STA[data_STA['flag']==4].resample('1min').mean() data_STA_5=data_STA[data_STA['flag']==5].resample('1min').mean() data_STA=pd.concat([data_STA_4,data_STA_5],axis=0) #合并两种标气 data_STA=data_STA.dropna() # data_sampleL = data[data['flag'] == 2].resample('5min').mean() data_STA.sort_index(inplace=True) data_STA.to_csv('data_STA.csv') print(data_STA) '''4.线性订证，输出高低层订证后的数据''' data = CONC(data_sampleH, data_STA) #日平均、月平均平均 data_sampleH_D=data.resample('D').mean() data_sampleH_D.to_csv('DZ_data_H_Day.csv') data_sampleH_M=data_sampleH_D.resample('M').mean() data_sampleH_M.to_csv('DZ_data_H_Month.csv')

订证后的分析

#Author Wu Dongqiao 2020.12.8 import pandas as pd import numpy as np import matplotlib.pyplot as plt #日数据 data=pd.read_csv(r'D:\data\picarro\DZ_data_H_Day.csv',index_col='time') data.index=pd.to_datetime(data.index) #月数据 df=pd.read_csv(r'D:\data\picarro\DZ_data_H_Month.csv',index_col='time') df.index=pd.to_datetime(df.index) # 作图段 fig1, axes = plt.subplots(nrows=3, figsize=(12, 10), sharex=True) # 设置一次画3个图 #将日与月画在一张图上 data['CO2'].plot(ax=axes[0], grid=True, legend=True) df['CO2'].plot(ax=axes[0], grid=True, legend=True) data['CO'].plot(ax=axes[1], grid=True, legend=True) df['CO'].plot(ax=axes[1], grid=True, legend=True) data['CH4'].plot(ax=axes[2], grid=True, legend=True) df['CH4'].plot(ax=axes[2], grid=True, legend=True) axes[0].set_ylabel('ppm') axes[1].set_ylabel('ppm') axes[2].set_ylabel('ppm') plt.tight_layout() plt.savefig('-温室气体变化.jpg') plt.show()