利用MATLAB求臭氧MDA8数据第90百分位数

文摘   其他   2024-07-23 16:26   广东  

利用MATLAB求

臭氧MDA8数据第90百分位数

——文末附完整代码


作者:第八星系-李智、成信大本科22级-陈嘉慧

邮箱:lizhi258147369@163.com



背景介绍

本文所用公式参考《环境空气质量评价技术规范HJ663-2013》附录A(规范性附录)数据统计方法,具体公式如下。






单年臭氧浓度的第90百分位数

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)为污染物浓度序列中的浓度值数量

end


臭氧浓度第90百分位数

kk=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是我们要的超标天数')

在命令行窗口中显示这句话来描述变量






逐月臭氧浓度的第90百分位数

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


臭氧浓度第90百分位数

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






工作区(以逐月为例)






命令行窗口(以逐月为例)



完整代码(单年)

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)为污染物浓度序列中的浓度值数量
end
kk=round(k);  % 将k的值四舍五入for i=1:nAenddata(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”

或扫码加入我们~


本文编辑:周游

第八星系人造大气理论爱好者
记录与交流python、matlab等科研工具。记录与交流大气科学的学科知识
 最新文章