MATLAB实现记录-CFANZ编程社区

1、太阳光度计&mod04

1.1光度计数据（day平均）相应的modis数据提取

剔除双方残缺值
整合数据到相应天数栏

clc;
clear;

data=xlsread('D:\study\AOD\2020beijing\data.xls')
P =[1,2,3,4];
j = 2;
Files = dir(fullfile('D:\study\AOD\2020beijing\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\2020beijing\data.xls', P);   %写入表格
p=xlsread('D:\study\AOD\2020beijing\data.xls');

1.2光度计数据相关性分析

data:不同波段影像数据
相关系数
regress()拟合后的数据残差图

clc;
clear;

data=xlsread('D:\study\AOD\2020beijing\data.xls');
data2=xlsread('D:\study\AOD\2020beijing\data2.xlsx');
data3=xlsread('D:\study\AOD\2020beijing\data3.xlsx');
y=data2(:,1);
y=y';
x=data2(:,2)
x=x';
B=data(:,3);
B=B';
b=data(:,6);
b=b';

blue=corrcoef(B,b); %相关系数
an=corrcoef(y,x);

% plot(x,y,'k*',x1,y1,'b--')

X = [ones(length(y),1), x'];%x'表示行向量转置为列向量
Y = y';
[ b,bint,r,rint,stats ] = regress(Y,X);
rcoplot(r,rint)
a=b(1)+(b(2))*x;
% plot(x,Y,'*',x,a,'-')

B=data(:,3);
B=B';
b=data(:,6);
b=b';
X = [ones(length(b),1), B'];%x'表示行向量转置为列向量
Y = b';
[ b,bint,r,rint,stats ] = regress(Y,X);
rcoplot(r,rint)
a=b(1)+(b(2))*B
plot(B,Y,'*',B,a,'-')

1.3 modis时间相应小时的太阳光度计数据提取

提取与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;%同一天数的时间的个数
n = 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)
                k = j;
                n = n+1;
                A = A+data(j,7);%level1和
                B = B+data(j,8);
            end
        end

end

s = 60;
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)
            k=j;
            n=n+1;
            A = A+data(j,7);
            B = B+data(j,8);   
        end
    end
end
end
J= J+a
k%选择数据的行数
A=A/n
B=B/n
D =horzcat(day,A,B,k);
P=vertcat(P,D);
dday=day
end

xlswrite ('D:\study\AOD\2021beijing\data6_6.xls', P);

1.5 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;%同一天数的时间的个数
n = 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)
                k = j;
                s = S
                A = data(j,7);%level1和
                B = data(j,8);
            end
        end

end

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)
            k=j;
            s = S
            A = data(j,7);
            B = data(j,8);   
        end
    end
end
end
J= J+a
k%选择数据的行数
A = data(K,7);
B = data(K,8); 
jj=K-3;
jjj=K+3
AA = mean(data(jj:jjj,7))
BB = mean(data(jj:jjj,8))
D =horzcat(day,AA,BB,k);
P=vertcat(P,D);
dday=day
end

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

1.6 写入有效的modis&&太阳光度计数据并给出相关系数

clc;
clear;

data1=xlsread('D:\study\AOD\2021beijing\data1.xls');
data2 = xlsread('D:\study\AOD\2021beijing\data2.xls');

P=[];
blie=[];
day=data2(:,1);
A=data2(:,2);%暗
B=data2(:,3);%蓝
C=data2(:,4);%结
for i=0:7
    n=6+i*3;
    m=7+i*3;
    nn=1+i*4;
    mm=2+i*4;
    mmm=3+i*4;
E = data2(:,n);
F = data2(:,m);

a=0
for j=2:70
    d=day(j);
    b=B(j);
    e=E(j);
    if ((~isnan(b)) && (~isnan(e)) && (e>0))
%         D= horzcat(d,b,e);
a=a+1
        P(a,nn)=d;
        P(a,mm)=b;
        P(a,mmm)=e;
    end
end

r=corrcoef(P(:,mm),P(:,mmm)) 
a=i+1
bule(a)=r(1,2)
end

xlswrite ('D:\study\AOD\2021beijing\data_blue.xls', P);