2019级 505组 源代码
最终的计算结果部分由python程序给出,另一部分由matlab计算得到
https://www.yuque.com/books/share/a0b1ee7d-cb7b-4b75-9bee-f8782e7b9ba8?
python程序如下
- 实验工具:
使用Python语言完成,主要应用numpy库和math库完成计算和存储 - 实验代码:
import numpy as np
import math
import pandas as pd
import os
path = r'D:\学习\大二\大探\涡度\05\25m'
llist = os.listdir(path)
filelist = []
for text in llist:
if text.endswith('.dat') == True:
filelist.append(text)
def wind_data_clean(start, end, n1, n2, n3):
U, V, W = [], [], []
u0 = np.array(df[n1][start:end])
v0 = np.array(df[n2][start:end])
w0 = np.array(df[n3][start:end])
l1 = np.nanmean(u0) + 3 * np.nanstd(u0)
r1 = np.nanmean(u0) - 3 * np.nanstd(u0)
U = df[n1][start:end]
l2 = np.nanmean(v0) + 3 * np.nanstd(v0)
r2 = np.nanmean(v0) - 3 * np.nanstd(v0)
V = df[n2][start:end]
l3 = np.nanmean(w0) + 3 * np.nanstd(w0)
r3 = np.nanmean(w0) - 3 * np.nanstd(w0)
W = df[n3][start:end]
for i in range(start , end):
if U[i] > l1 or U[i] < r1:
U[i] = np.nan
V[i] = np.nan
W[i] = np.nan
for i in range(start , end):
if V[i] > l2 or V[i] < r2:
U[i] = np.nan
V[i] = np.nan
W[i] = np.nan
for i in range(start , end):
if W[i] > l3 or W[i] < r3:
U[i] = np.nan
V[i] = np.nan
W[i] = np.nan
U = np.array(U)
V = np.array(V)
W = np.array(W)
return U, V, W
def change_NAN(start, end, name):
for i in range(start, end):
if df[name][i] == 999:
df[name][i] = np.nan
return 0
def data_clean(start, end, name):
u0 = np.array(df[name][start:end])
l = np.nanmean(u0) + 3 * np.nanstd(u0)
r = np.nanmean(u0) - 3 * np.nanstd(u0)
for i in range(start , end):
if df[name][i] > l or df[name][i] < r:
df[name][i] = np.nan
U = 0
return U
def average_data(start, end, data_):
average_data_ = np.nanmean(data_[start:end])
return average_data_
def change_xyz(start, end, Ux_, Uy_, Uz_):
u0 = Ux_[start:end]
v0 = Uy_[start:end]
w0 = Uz_[start:end]
cos_alpha = np.nanmean(u0)/math.sqrt(np.nanmean(u0) ** 2 + np.nanmean(v0) ** 2)
sin_alpha = np.nanmean(v0)/math.sqrt(np.nanmean(u0) ** 2 + np.nanmean(v0) ** 2)
u1 = u0*cos_alpha + v0 * sin_alpha
v1 = -u0*sin_alpha + v0 * cos_alpha
w1 = w0
cos_beta = np.nanmean(u1)/math.sqrt(np.nanmean(u1) ** 2 + np.nanmean(w1) ** 2)
sin_beta = np.nanmean(w1)/math.sqrt(np.nanmean(u1) ** 2 + np.nanmean(w1) ** 2)
u2 = u1*cos_beta + w1 * sin_beta
v2 = v1
w2 = -u1*sin_beta + w1 * cos_beta
return u2, v2, w2
def ganre(start, end):
p0 = df['amb_press'][start:end]*1000
t0 = df['Ts'][start:end]+273.15
H2O = df['H2O'][start:end]/1000
rou0 = p0/(287.059*t0) + H2O
cp = 1004.67*(1+0.84*H2O/rou0)
t0 = np.array(t0)
Uz_ = np.array(Uz)
WT = np.nanmean(np.multiply((t0 - np.nanmean(t0)),(Uz_ - np.nanmean(Uz_))))
H = rou0*cp*WT
H = np.array(H)
return np.mean(H)
def qianre(start, end):
t0 = df['Ts'][start:end]+273.15
lamda = (2.501-0.00237*t0)*(10**6)
H2O = df['H2O'][start:end]/1000
Uz_ = np.array(Uz)
H2O_ = np.array(H2O)
WR = np.nanmean((H2O_ - np.nanmean(H2O_))*(Uz_ - np.nanmean(Uz_)))
return np.mean(lamda*WR)
def wind_direction(start, end):
u0 = Ux
v0 = Uy
cos_alpha = np.nanmean(u0)/math.sqrt(np.nanmean(u0) ** 2 + np.nanmean(v0) ** 2)
direction = math.acos(cos_alpha)*180/math.pi
return direction
def CO(start, end):
CO2 = df['CO2'][start:end]*0.000001
Uz_ = np.array(Uz)
CO2_ = np.array(CO2)
WC = np.nanmean((CO2_ - np.nanmean(CO2_))*(Uz_ - np.nanmean(Uz_)))
return np.mean(WC)
data = []
H_data_item = []
E_data_item = []
direction_item = []
C_data_item = []
for x in range(0,len(filelist)):
df = pd.read_csv(filelist[x],
header=None,
names=['TIMESTAMP','RECORD','Ux','Uy','Uz','Ts','diag_sonic','CO2','H2O','diag_irga','amb_tmpr','amb_press','CO2_sig_strgth','H2O_sig_strgth'],
skiprows=lambda x: (x%36000)//600<30,
engine='python')
aUx, aUy, aUz = [], [], []
H_data = []
E_data = []
wind_d = []
C_data = []
for i in range(0, len(df["Ux"]), 18000):
a = change_NAN(i, i+18000, 'Ux')
a = change_NAN(i, i+18000, 'Uy')
a = change_NAN(i, i+18000, 'Uz')
a = change_NAN(i, i+18000, 'Ts')
a = change_NAN(i, i+18000, 'H2O')
a = change_NAN(i, i+18000, 'amb_press')
a = data_clean(i, i+18000, 'CO2')
Ux, Uy, Uz = wind_data_clean(i, i+18000, 'Ux', 'Uy', 'Uz')
wind_d.append(wind_direction(i, i+18000))
Ux, Uy, Uz = change_xyz(0, 18000, Ux, Uy, Uz)
aUx.append(average_data(0, 18000, Ux))
aUy.append(average_data(0, 18000, Uy))
aUz.append(average_data(0, 18000, Uz))
H_data.append(ganre(i, i+18000))
E_data.append(qianre(i, i+18000))
C_data.append(CO(i, i+18000))
H_data_item.append(H_data)
E_data_item.append(E_data)
direction_item.append(wind_d)
C_data_item.append(C_data)
data.append(H_data_item)
data.append(E_data_item)
data.append(direction_item)
data.append(C_data_item)
np.save('25m_10_13-15.npy', data)
matlab相关程序
1.DealData.m
clc;clear;close all;
Len=18000; % 每30min数据有18000行
data0temp=importdata("D:\exp\exp3\05\data\25m1013.dat") ; % 加载数据 importdata
data0=data0temp.data; % data0为存放原始数据的矩阵
for ii=1:48 %1day - 24hours - 48 halfhours
Rawdata0(:,1:4)=data0(1+(ii-1)*Len:Len*ii,2:5);
% 2,3,4列为风速在3个方向上的分量,第5列为气温数据 在Rawdata数据中,1-3列存风速数据
Rawdata0(:,5:6)=data0(1+(ii-1)*Len:Len*ii,7:8);
% 第7列数据为CO2密度,第8列为水汽密度D
%进行对数据的质量控制,依旧3倍标准差原则筛选*******************************************************
clear ave; clear std_data;
for j=1:6
ave(j) = mean(Rawdata0(:,j),'omitnan');
std_data(j) = std(Rawdata0(:,j),'omitnan');
end
ave = ave'; std_data = std_data';
for j=1:6
for i= 1:18000
if Rawdata0(i,j)>ave(j)+3*std_data(j)|Rawdata0(i,j)<ave(j)-3*std_data(j)
Rawdata0(i,j) = nan;
end
end
end
clear Rawdata_1
Logi2=isnan(Rawdata0(:,1))|isnan(Rawdata0(:,2))|isnan(Rawdata0(:,3))|isnan(Rawdata0(:,4))|isnan(Rawdata0(:,5))|isnan(Rawdata0(:,6));
for j=1:6
Rawdata_1(:,j)=Rawdata0(~Logi2,j); % 剔除数据中含有NaN的行(该行对应的布尔值为0)
end
density_of_air = 1.293;
%进行坐标旋转*******************************************************************************************
data_1=rotate1(Rawdata_1);
datamean0_1(ii,:)=mean(Rawdata_1);
% % mean函数作用于矩阵求的就是各列的平均值
datamean_1(ii,:)=mean(data_1);
% 计算旋转后平均值,以datamean_1的结果为准!!!
%计算协方差矩阵*******************************************************************************
datacov_1=zeros(2,0);
datacov_1=cov(data_1); % data_1为坐标旋转后的最终矩阵
%datacov 为协方差方针6*6
flux_1(ii,:)=datacov_1(3,:); % 垂直方向上的通量
sigma_1(ii,1)=datacov_1(1,1)^0.5; sigma_1(ii,2)=datacov_1(2,2)^0.5;
sigma_1(ii,3)=datacov_1(3,3)^0.5; sigma_1(ii,4)=datacov_1(4,4)^0.5;
% sigma_1存储标准差数据
%计算湍流强度,摩擦速度,动量通量**********************************************************
U(ii)=sqrt( datamean_1(ii,1)^2+datamean_1(ii,2)^2+datamean_1(ii,3)^2 );
%calculate average wind speed
for j=1:3
tur_inten(ii,j) = sigma_1(ii,j)/U(ii);
end
U = U';
%计算磨擦速度 用fric_speed表示
fric_speed(ii) = (datacov_1(1,3)^2+datacov_1(2,3)^2)^0.25;
% 计算摩擦速度,也可直接利用垂直方向通量数据:fric_speed(ii) = (flux_1(ii,1)^2+flux_1(ii,2)^2)^0.25;
%计算切应力(动量通量),用mom_flux表示
mom_flux = -density_of_air*fric_speed;
end
fric_speed = fric_speed';
mom_flux = mom_flux';
mean_wind = datamean_1(2:2:48,1:3);
for i = 1:24
frictionspeed(i) = fric_speed(2*i);
tur(i,:) = tur_inten(2*i,:);
windstd(i,1:3) = sigma_1(2*i,1:3);
momen(i) = mom_flux (2*i);
end
frictionspeed = frictionspeed';
momen = momen';
flux_1=flux_1(2:2:48,:);
2.rotate1.m
function data=rotate1(Rawdata)
%该程序的目的是做坐标变换,TR
%input:Rawdata
%output:data
%
% 只通过旋转改变传入参数矩阵前3列的风速数据,不影响第4-6列的其他数据
datamean=mean(Rawdata);
U=datamean(1);V=datamean(2);W=datamean(3);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%旋转,使得U沿着平均速度方向,而V=0,
YawMatrix=[U/(sqrt(U^2+V^2)),-V/(sqrt(U^2+V^2)),0;V/(sqrt(U^2+V^2)),U/(sqrt(U^2+V^2)),0;0,0,1];
YawData=Rawdata(:,1:3)*YawMatrix;
datamean=mean(YawData);
U=datamean(1);V=datamean(2);W=datamean(3);
%通过乘变换方阵实现旋转
%旋转,使得U沿着平均速度方向,而W=0,
PitchMatrix=[U/(sqrt(U^2+W^2)),0,-W/(sqrt(U^2+W^2)); 0,1,0;W/(sqrt(U^2+W^2)),0,U/(sqrt(U^2+W^2))];
PitchData=YawData*PitchMatrix;
datamean=zeros(2,0);datamean=mean(PitchData);
datacov=zeros(2,0);datacov=cov(PitchData);
data=Rawdata;
data(:,1:3)=PitchData;
3.程序A(基于向量的数据处理程序)
DealDataFormer.m
i=1;
for(j=18001:36000:846001)
Ux_raw_ave(i) = mean(Ux_raw(j:j+17999),'omitnan');
Ux_raw_std(i) = std(Ux_raw(j:j+17999),'omitnan');
Uy_raw_ave(i) = mean(Uy_raw(j:j+17999),'omitnan');
Uy_raw_std(i) = std(Uy_raw(j:j+17999),'omitnan');
Uz_raw_ave(i) = mean(Uz_raw(j:j+17999),'omitnan');
Uz_raw_std(i) = std(Uz_raw(j:j+17999),'omitnan');
i=i+1;
end
Ux = Ux_raw;
Uy = Uy_raw;
Uz = Uz_raw;
% 后续计算只需要每小时后30分钟的数据,故只对这30分钟数据进行筛选和修正,因此循环不遍历全部数据(如下)
for i=1:24
for j=36000*i-17999:36000*i
if (Ux_raw(j)>Ux_raw_ave(i)+3*Ux_raw_std(i) || Ux_raw(j)<Ux_raw_ave(i)-3*Ux_raw_std(i))
Ux(j)=NaN;
Uy(j)=nan; Uz(j)=nan;
end
if (Uy_raw(j)>Uy_raw_ave(i)+3*Uy_raw_std(i) || Uy_raw(j)<Uy_raw_ave(i)-3*Uy_raw_std(i))
Uy(j)=NaN;
Ux(j) = nan; Uz(j) = nan;
end
if (Uz_raw(j)>Uz_raw_ave(i)+3*Uz_raw_std(i) || Uz_raw(j)<Uz_raw_ave(i)-3*Uz_raw_std(i))
Uz(j)=NaN;
Ux(j)=nan; Uy(j)=nan;
end
end
end
for i=1:24
for j=36000*i-17999:36000*i
if isnan(Ux(j)) | isnan(Uy(j)) |isnan(Uz(j))
Ux(j) = nan;
Uy(j) = nan;
Uz(j) = nan;
end
end
end
%重新计算平均值,生成向量
i = 1;
for(j=18001:36000:846001)
Ux_ave(i) = mean(Ux(j:j+17999),'omitnan');
Ux_std(i) = std(Ux(j:j+17999),'omitnan');
Uy_ave(i) = mean(Uy(j:j+17999),'omitnan');
Uy_std(i) = std(Uy(j:j+17999),'omitnan');
Uz_ave(i)= mean(Uz(j:j+17999),'omitnan');
Uz_std(i) = std(Uz(j:j+17999),'omitnan');
i=i+1;
end
% 进行第1次坐标旋转
%
Ux1 = Ux; Uy1 = Uy;
i = 1;
while(i<=24)
alpha(i)= atan(Uy_ave(i)/Ux_ave(i));
Ux1(i*36000-17999:i*36000) = Ux(i*36000-17999:i*36000) * cos(alpha(i))+ Uy(i*36000-17999:i*36000)*sin(alpha(i));
Uy1(i*36000-17999:i*36000) =-1*Ux(i*36000-17999:i*36000) * sin(alpha(i))+ Uy(i*36000-17999:i*36000)*cos(alpha(i));
i = i+1;
end
%
% 进行第2次坐标旋转
%
% 先算第一次坐标旋转后Ux的平均值(24个时间段,每段半小时)
i=1;
while(i<=24)
Ux1_ave(i) = mean(Ux1(i*36000-17999:i*36000),'omitnan');
i=i+1;
end
Ux2 = Ux1; Uz2 = Uz;
i =1 ;
while(i<=24)
beta(i) = atan(Uz_ave(i)/Ux1_ave(i));
Ux2(i*36000-17999:i*36000) = Ux1(i*36000-17999:i*36000)*cos(beta(i))+Uz(i*36000-17999:i*36000)*cos(beta(i));
Uz2(i*36000-17999:i*36000) = -Ux1(i*36000-17999:i*36000)*sin(beta(i))+Uz(i*36000-17999:i*36000)*cos(beta(i));
i=i+1;
end
% 计算坐标旋转后风速的平均值,作为最终结果
i = 1;
while(i<=24)
Ux_ave_final(i) = mean(Ux2(i*36000-17999:i*36000),'omitnan');
Uy_ave_final(i) = mean(Uy1(i*36000-17999:i*36000),'omitnan');
Uz_ave_final(i) = mean(Uz2(i*36000-17999:i*36000),'omitnan');
Ux_std_final(i) = std(Ux2(i*36000-17999:i*36000),'omitnan');
Uy_std_final(i) = std(Uy1(i*36000-17999:i*36000),'omitnan');
Uz_std_final(i) = std(Uz2(i*36000-17999:i*36000),'omitnan');
i=i+1;
end
i = 1;
while(i<=24)
uv(i) = mean((Ux2(i*36000-17999:i*36000)-Ux_ave_final(i)).*(Uz2(i*36000-17999:i*36000)-Uz_ave_final(i)),'omitnan');
i=i+1;
end
U_average = [Ux_ave_final;Uy_ave_final;Uz_ave_final]';
U_std = [Ux_std_final;Uy_std_final;Uz_std_final]';
%计算风向
i = 1;
while(i<=24)
judge = Uy_ave_final(i)/Ux_ave_final(i);
if (Uy_ave_final(i)>0)
D(i) = 270 - atand(judge);
end
if (Uy_ave_final(i)<0)
D(i) = 90 - atand(judge);
end
if ((Ux_ave_final(i)>0) && (Uy_ave_final(i)==0))
D(i) = 180;
end
if ((Ux_ave_final(i)<0) && (Uy_ave_final(i)==0))
D(i) = 0;
end
i=i+1;
end