利用MATLAB求
臭氧MDA8数据第90百分位数
——文末附完整代码
作者:第八星系-李智、成信大本科22级-陈嘉慧
邮箱:lizhi258147369@163.com
完整代码(单年)
clear
clc
cdata1 = xlsread('cities_O3.xlsx','2022','B2:AO367');
% 从工作簿指定工作表“2022”的指定范围内读取数据
datasize=size(cdata1);
% 返回以元素是cdata1的相应维度的长度的行向量
h=datasize(1);
n=datasize(2);
[data2,index]=sort(cdata1);
% 根据下图要求,将污染物浓度序列按数值由小到大排序
for m=1:n
a=data2(:,m);
NaNnumbel(m)=numel(find(isnan(a)));
NotNaNnumbel(m)=h-NaNnumbel(m);
k(m)=1+(NotNaNnumbel(m)-1)*0.9;
% 根据公式,计算90%对应的序数k
% NotNaNnumbel(m)为污染物浓度序列中的浓度值数量
endkk=round(k); % 将k的值四舍五入
for i=1:n
Aenddata(i)=data2(kk(i),i);
end
% 找序数对应的浓度
disp('工作区的‘Aendddate’是我们要的臭氧浓度90百分位数')
% 在命令行窗口中显示这句话来描述变量
for i=1:n
nongdu=cdata1(:,i);
Aasum(i)=sum(sum(nongdu>160));
end
disp('工作区的‘Aasum’是我们要的超标天数')
% 在命令行窗口中显示这句话来描述变量
1
完整代码(逐月)
clear
clc
cdata1 = xlsread('cities_O3.xlsx','2022','B2:AP366');
% 从工作簿指定工作表“2022”的指定范围内读取数据
datasize = size(cdata1);
% 返回以元素是cdata1的相应维度的长度的行向量
h = datasize(1);
n = datasize(2);
cdata2 = mat2cell(cdata1, [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]);
% 将数组 cadata1 划分为一个 a×1 元胞数组,其中 a 等于方括号中元素的数量
%将数组 cadata1 划分为更小的数组,并在元胞数组 cadata2 中返回它们
months = {'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'};
Aenddata = zeros(12, n);
Aasum = zeros(12, n);
% 返回一个 12×n 矩阵
for month = 1:12
[month_data,~] = sort(cdata2{month});
datasize = size(month_data);
h = datasize(1);
% 重复1-12月的for循环
% 根据公式计算序数
NaNnumbel = zeros(1, n);
NotNaNnumbel = zeros(1, n);
k = zeros(1, n);
% 返回一个 1×n 矩阵
for m = 1:n
a = month_data(:, m);
NaNnumbel(m) = numel(find(isnan(a)));
NotNaNnumbel(m) = h - NaNnumbel(m);
k(m) = 1 + (NotNaNnumbel(m) - 1) * 0.9;
% 套用公式,与单年计算同理
end
kk = round(k); % 四舍五入
for i = 1:n
Aenddata(month, i) = month_data(kk(i), i);
end
disp(['工作区的''Aendddate''即臭氧浓度90百分位数 - ' months{month}]);
% 在命令行窗口中显示这句话来描述变量
for i = 1:n
nongdu = month_data(:, i);
Aasum(month, i) = sum(sum(nongdu > 160));
end
disp(['工作区的''Aasum''即超标天数 - ' months{month}]);
% 在命令行窗口中显示这句话来描述变量
end
第八星系
如需获取数据,
请后台回复“臭氧MDA8”
或扫码加入我们~
本文编辑:周游