太阳光度计CE-318数据处理

太阳光度计CE-318数据处理

备注：处理公式

在我国近海，α的值在0到3之间，所以他们相对误差最大不超过25%，而通过查阅相关资料，北京地区α的值可以近似的取1.665。

大气是不断运动的，气溶胶在短时间内也可能会发生很大的变化，而MODIS传感器和AERONET观测站并不能保证在同一时刻获取数据，所以，经过查阅相关资料，根据Ichoku等人的研究结论，MODIS数据从30km×30km到90km×90km范围内变化对窗口平均值的影响比较小[52]，而Zhao等人也认为100km/上1h的时空匹配窗口是最合适的[53。因此，本文对AERONET观测站周围90km×90km的区域进行空间采样，取MODIS卫星遥感数据前后1h范围内的AERONET观测数据的平均值进行反演结果的对比。

一、数据下载

下载链接：https://aeronet.gsfc.nasa.gov/cgi-bin/draw_map_display_aod_v3

在右侧地图中选择需要的站点、等级数据

二、数据处理

2.1太阳光度计550nm数据计算

y:550nm
x:500nm
y=[(550/500)^(-ɑ)]* x

2.2太阳光度计时间转换

[~,~,raw]=xlsread('D:\study\AOD\2021beijing\550.xls');
%[num,txt,raw]num数值，txt字符串，raw数值与字符串
data=raw(1:10900,:);

时间hh:mm:ss改为（自定义选项）hhmm

提取的hhmm数字转文本

字符串提取hh和mm：

 =TEXT(LEFT(F3,2),"00")     前两个字符
 =TEXT(RIGHT(F3,2),"00")    后两个字符
 =TEXT(MIN(F3,3,2),"00")  %%中间 
 转为数值

a指数

=-((LN(D2/F2))/(LN(440/870)))
=((550/500)^(-G2))*E2   AOD_550nm_500
=((550/440)^(-G2))*D2   AOD_550nm_440
=((550/500)^(-1.665))*E2    AOD_550nm_固定

（2）数据处理代码.m

clc;
clear;
%%%将太阳光度计站点原始数据提取，加工
%%%1、day 2、h,m的转换 3、440,500,870nm数据 4、a指数 5、550nm数据

%%%更改路径csv_path即可  D:\study\CE-318\北京N39.977,E116.381\
csv_path=  'D:\study\CE-318\北京-CAMS(39.933牛顿，116.317牛顿)\';   %文件夹路径 
path_list = dir(strcat(csv_path,'*.csv'));  %dir 函数 列出当前目录下所有子文件夹和文件%
list_num = length(path_list);%%文件数量
Day = 1;  %天数
P = 2;    %小时
p = 3;    %分钟
A = 4;    %870nmAOD
B = 5;    %500nmAOD
C = 6;    %440nmAOD
D = 7;    %a指数
E = 8;    %500nm的550nmAOD
F = 9;    %440nm的550nmAOD
G = 10; %北京固定a指数550nmAOD

for i=1:list_num
    filename = [csv_path,path_list(i).name];  %获取文件信息，以便最后取名
    yearname = path_list(i).name(1:4);
    AAA = [];   %每一次循环清空AAA
    delete('AAA.xls'); 
    AAA = take(filename,[1, Inf]);   %take为函数，提取太阳光度计时间，数值等原始信息
    writetable(AAA, 'AAA.xls');      %提取的AAA数据为table数据，writetable写入为xls
    [~,~,raw] = xlsread('AAA.xls');  %读取xls数据，包括字符串
    raw(2:8,:)=[];    %删除前几行空白      
    raw(1,:)={0,0,0,0,0};  %取表格头为0 
    time= raw(:,1);        
    day = raw(:,2);        %%%%分别获取相应列数据
    a = raw(:,3);
    b = raw(:,4);
    c = raw(:,5);
    
    [K,n]=size(time);  %获取数据量K
    for k=2:K   %循环每一行数据
    T = time{k};%%%读取cell时{}为数值；（）为cell原数据类型
    dd = day{k};
    aa = a{k};
    bb = b{k};
    cc = c{k};
    
    if length(T)>17    %字符串长度length
    H = T(11:12);      %根据字符串长度判断，获取时间信息，H小时，M分钟
    M = T(14:15);
    H = str2num(H);
    M = str2num(M);
    elseif length(T)>9
        H = T(11);
        M = T(13:14);
        H = str2num(H);
        M = str2num(M);
    else
        H = 0;
        M = 0;
    end
    P = vertcat(P,H);        %将每一次数据垂直方向拼接
    p = vertcat(p,M);
    
    DD = -((log10(cc/aa))/(log10(440/870)));
    EE = ((550/500)^(-DD))*bb;%%500的
    FF = ((550/440)^(-DD))*cc;%%440的
    GG = ((550/500)^(-1.665))*bb;%%固定
    
    Day = vertcat(Day,dd);   %将每一次数据垂直方向拼接
    A = vertcat(A,aa);
    B = vertcat(B,bb);
    C = vertcat(C,cc);
    D = vertcat(D,DD);
    E = vertcat(E,EE);
    F = vertcat(F,FF);
    G = vertcat(G,GG);
    end
    DATA = horzcat(Day,P,p,A,B,C,D,E,F,G); %水平拼接
%     DATA = horzcat(Day,P,p,A,B,C,D,E,F);
    Day = 1;
    P = 2;
    p = 3;
    A = 4;
    B = 5;
    C = 6;
    D = 7;
    E = 8;
    F = 9;
    G = 10;
    xlswrite(strcat(csv_path,yearname,'.xlsx'),DATA); %自定义命名将DATA写入表格
end

2.3 制作表格：太阳光度计数据天数时间对标modis天数时间

clc;
clear;
%制作表格：太阳光度计数据天数时间对标modis天数时间

data = xlsread('D:\study\AOD\2021beijing\data.xls');
data = data(:,1:10);
sunday = data(:,2);
sunhour =data(:,9);
sunss = data(:,10);
P =[1,2,3,4];
J=1;
a=0

Files = dir(fullfile('D:\study\AOD\2021beijing\tiff\*.tif')); % 读取文件夹内的mat格式的文件
LengthFiles = length(Files); %所有文件的数量


% for i=1:LengthFiles

for i=1:70 
name=Files(i).name;           %读取struct变量的格式
day = name(15:17);
day = str2num(day);   %字符串转数字
hour = name(19:20);
hour = str2num(hour);
ss = name(21:22);
ss  = str2num(ss);

A = 0;
B = 0;

if(sunday(J)~=day)
K = J-a
I = i+1
Name=Files(I).name;           %读取struct变量的格式
Day = Name(15:17);
Day = str2num(Day);   %字符串转数字
if (i>1 && dday==day)
    J=J-a;
elseif(sunday(J)==Day)
    k=-1;
    D =horzcat(day,A,B,0);
    P=vertcat(P,D);
    continue
end
end

    
s = 60;
S = 0;
a = 0;
for j=J:J+100
    if(sunday(j)==day)
        a=a+1;
    end
end
for j=J:J+a-1
        if(sunhour(j)==hour)
            S = abs(sunss(j)-ss);
            if(S<s)
                s = S;
                k = j;
                A = data(j,7);%level1
                B = data(j,8);
            end
        end
end
J
A
a

s = 240;
S = 0;
if (A==0)
for j=J:J+a-1
    if(sunhour(j)~=hour)
        H = abs(sunhour(j)-hour)-1
        if (sunhour(j)<hour)
        S=H*60+(60-sunss(j)+ss);
        else
            S= H*60+(60-ss+sunss(j));
        end
        if(S<s)
            s = S;
            k=j;
            A = data(j,7);
            B = data(j,8);   
        end
    end
end
end
J= J+a
k%选择数据的行数
D =horzcat(day,A,B,k);
P=vertcat(P,D);
dday=day
end

xlswrite ('D:\study\AOD\2021beijing\data1.xls', P);
p=xlsread('D:\study\AOD\2021beijing\data1.xls');

2.4 提取modis数据对应到太阳光度计数据

clc;
clear;

data=xlsread('D:\study\AOD\2021beijing\data1.xls');
P =[1,2,3,4];
j = 2;
Files = dir(fullfile('D:\study\AOD\2021beijing\TIFF\*.tif')); % 读取文件夹内的mat格式的文件
LengthFiles = length(Files); %所有文件的数量

for i=1:LengthFiles
name=Files(i).name;           %读取struct变量的格式
day = name(15:17);
day = str2num(day);   %字符串转数字
folder=Files(i).folder;
[Data,R] = geotiffread(strcat(folder,'\',name)); 
    
% [Data,R] = geotiffread('D:\study\AOD\158.tif');

row = R.RasterSize(1);
col = R.RasterSize(2);
X = (R.LatitudeLimits(2)-39.977)/R.RasterExtentInLatitude;
Y = (116.381-R.LongitudeLimits(1))/R.RasterExtentInLongitude;
x = floor(X*row);  %指定经纬度所在行列
y = floor(Y*col);

a = Data(:,:,1);
A = a(x,y);
b = Data(:,:,2);
B = b(x,y);
c = Data(:,:,3);
C = c(x,y);
D = horzcat(day,A,B,C);
if (D(1)==data(j-1))
    P=vertcat(P,D);
end
if (D(1)==data(j))
P=vertcat(P,D);
j=j+1;
end
end

xlswrite ('D:\study\AOD\2021beijing\data2.xls', P);
p=xlsread('D:\study\AOD\2021beijing\data2.xls');