太阳光度计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指数550nmAODfor i=1:list_numfilename = [csv_path,path_list(i).name]; %获取文件信息，以便最后取名yearname = path_list(i).name(1:4);AAA = []; %每一次循环清空AAAdelete('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); %获取数据量Kfor k=2:K %循环每一行数据T = time{k};%%%读取cell时{}为数值；（）为cell原数据类型dd = day{k};aa = a{k};bb = b{k};cc = c{k};if length(T)>17 %字符串长度lengthH = T(11:12); %根据字符串长度判断，获取时间信息，H小时，M分钟M = T(14:15);H = str2num(H);M = str2num(M);elseif length(T)>9H = T(11);M = T(13:14);H = str2num(H);M = str2num(M);elseH = 0;M = 0;endP = 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);endDATA = 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=0Files = dir(fullfile('D:\study\AOD\2021beijing\tiff\*.tif')); % 读取文件夹内的mat格式的文件
LengthFiles = length(Files); %所有文件的数量% for i=1:LengthFilesfor 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
ends = 60;
S = 0;
a = 0;
for j=J:J+100if(sunday(j)==day)a=a+1;end
end
for j=J:J+a-1if(sunhour(j)==hour)S = abs(sunss(j)-ss);if(S<s)s = S;k = j;A = data(j,7);%level1B = data(j,8);endend
end
J
A
as = 240;
S = 0;
if (A==0)
for j=J:J+a-1if(sunhour(j)~=hour)H = abs(sunhour(j)-hour)-1if (sunhour(j)<hour)S=H*60+(60-sunss(j)+ss);elseS= H*60+(60-ss+sunss(j));endif(S<s)s = S;k=j;A = data(j,7);B = data(j,8); endend
end
end
J= J+a
k%选择数据的行数
D =horzcat(day,A,B,k);
P=vertcat(P,D);
dday=day
endxlswrite ('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
endxlswrite ('D:\study\AOD\2021beijing\data2.xls', P);
p=xlsread('D:\study\AOD\2021beijing\data2.xls');
