CALIPSO和AEOLUS卫星数据矩阵维度统一

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

本文用于记录分享处理CALIPSO和AEOLUS卫星数据的一些问题和解决。

一、矩阵维度统一

AEOLUS有24个高度层而CALIPSO有399个高度层，因此处理起来会有很多不便。所以该函数用于统一AEOLUS和CALIPSO的矩阵维度。
使用工具：MATLAB2022a
数据CALIPSO L2 aerosol profile & Aeolus L2A

二、1.使用步骤

使用前请更改CODA路径和文件根目录并选定经纬度
addpath ‘C:\Program Files\CODA\share\coda\definitions’
root_path_AEOLUS
root_path_CALIPSO
代码如下（示例）：

%%
%本文件用于调整AEOLUS和CALIPSO的矩阵的size以满足下一步需求
function [AEOLUS_altitude,AEOLUS_latitude,AEOLUS_longitude,...
   CALIPSO_altitude, CALIPSO_latitude,CALIPSO_longitude,...
   Backscatter_532,Backscatter_1064,Extinction_355,Extinction_532,...
   Extinction_1064] = resize(file_name_CALIPSO,file_name_AEOLUS)
%清理工作区&&CODA定义path
%clear all;clc; %#ok<*CLALL>
addpath 'C:\Program Files\CODA\share\coda\definitions'

%%
%文件信息
%根目录
root_path_AEOLUS=('D:\graduation_paper\graduation papar\data\Aeolus\in_area\');
root_path_CALIPSO=('D:\graduation_paper\graduation papar\data\CALIPSO\');

%文件名
% file_name_AEOLUS='AE_OPER_ALD_U_N_2A_20200914T225411022_005412007_011957_0004.DBL';
% file_name_CALIPSO='CAL_LID_L2_05kmAPro-Standard-V4-21.2020-09-14T18-39-46ZD.hdf';

file_path_AEOLUS=[root_path_AEOLUS,file_name_AEOLUS];
file_path_CALIPSO=[root_path_CALIPSO,file_name_CALIPSO];
%筛选范围
LonLim=[-150,-80];
LatLim=[20,60];

clear root_path_AEOLUS root_path_CALIPSO file_name_AEOLUS file_name_CALIPSO;

%%
%打开AEOLUS的消光系数和经纬度 
CODA_FILE_ID = coda_open(file_path_AEOLUS);
extinction_355_Origin=coda_fetch(CODA_FILE_ID,'sca_optical_properties',-1, ...
    'sca_optical_properties',-1,'extinction');

longitude_AEOLUS_Origin=coda_fetch(CODA_FILE_ID,'sca_optical_properties',-1, ...
    'geolocation_middle_bins',-1,'longitude');

latitude_AEOLUS_Origin=coda_fetch(CODA_FILE_ID,'sca_optical_properties',-1, ...
    'geolocation_middle_bins',-1,'latitude');

altitude_AEOLUS_Origin=coda_fetch(CODA_FILE_ID,'sca_optical_properties',-1, ...
    'geolocation_middle_bins',-1,'altitude');

%取AEOLUS经纬度的
AEOLUS_longitude=zeros(size(longitude_AEOLUS_Origin));
AEOLUS_latitude=zeros(size(latitude_AEOLUS_Origin));
for i=1:length(AEOLUS_latitude)
    temp=0;count=0;temp_Lon=cell2mat(longitude_AEOLUS_Origin(i));
    %先计算经度平均值
    for j=1:length(temp_Lon)
        if (temp_Lon(j)~=(-1.000000000000000e-06))
            temp=temp+temp_Lon(j);
            count=count+1;
        end
    end
    AEOLUS_longitude(i)=temp/count;
    
    temp=0;count=0;temp_Lat=cell2mat(latitude_AEOLUS_Origin(i));
    %后计算纬度平均值
    for j=1:length(temp_Lat)
        if (temp_Lat(j)~=(-1.000000000000000e-06))
            temp=temp+temp_Lat(j);
            count=count+1;
        end
    end
    AEOLUS_latitude(i)=temp/count;

end

%高度矩阵
AEOLUS_altitude=cell2mat(altitude_AEOLUS_Origin);
AEOLUS_altitude=reshape(AEOLUS_altitude,[24,length(altitude_AEOLUS_Origin)]);
AEOLUS_altitude=(AEOLUS_altitude');
AEOLUS_altitude(AEOLUS_altitude==-1)=NaN;
AEOLUS_altitude=AEOLUS_altitude./1000;
clear altitude_AEOLUS_Origin;

%矫正AEOLUS的经度计算方法(if>180 do-360)
 for i=1:length(AEOLUS_longitude)
     if(AEOLUS_longitude(i)>180.0)
         AEOLUS_longitude(i)= AEOLUS_longitude(i)-360.0;
     end
 end

clear longitude_AEOLUS_Origin latitude_AEOLUS_Origin 
clear temp count temp_Lat temp_Lon i j CODA_FILE_ID;


%整理AEOLUS消光系数
Extinction_355=cell2mat(extinction_355_Origin);
Extinction_355=reshape(Extinction_355,[24,length(extinction_355_Origin)]);
Extinction_355=Extinction_355';
clear extinction_355_Origin;

%%
%打开CALIPSO经纬度并计算平均值
CALIPSO_longitude=hdfread(file_path_CALIPSO,'Longitude');
CALIPSO_longitude=(CALIPSO_longitude(:,1)+CALIPSO_longitude(:,2)+CALIPSO_longitude(:,3))/3;
CALIPSO_latitude=hdfread(file_path_CALIPSO,'Latitude');
CALIPSO_latitude=(CALIPSO_latitude(:,1)+CALIPSO_latitude(:,2)+CALIPSO_latitude(:,3))/3;

%高度-0.5-20.2 60m/bin 20.2-30.1 180m/bin
Fixed_CALIPSO_altitude=-0.5:0.06:20.2;Fixed_CALIPSO_altitude(1)=[];
temp=20.2:0.18:30.1;temp(56)=[];temp(1)=[];Fixed_CALIPSO_altitude=[Fixed_CALIPSO_altitude,temp];clear temp;
Fixed_CALIPSO_altitude(:,:)=Fixed_CALIPSO_altitude(:,end:-1:1);


%读取后向散射系数消光系数
Backscatter_532=hdfread(file_path_CALIPSO,'Total_Backscatter_Coefficient_532');
Extinction_532=hdfread(file_path_CALIPSO,'Extinction_Coefficient_532');
Backscatter_1064=hdfread(file_path_CALIPSO,'Backscatter_Coefficient_1064');
Extinction_1064=hdfread(file_path_CALIPSO,'Extinction_Coefficient_1064');

%读取完毕清理工作区
clear file_path_CALIPSO file_path_AEOLUS


%%筛选经纬度，CALIPSO适当放宽半个经纬度防止算法溢出

%筛选AEOLUS经纬度在研究区域的坐标点
Origin_Length=length(AEOLUS_longitude);
for i=1:Origin_Length
    j=Origin_Length-i+1;
    if(AEOLUS_longitude(j)<LonLim(1)||AEOLUS_longitude(j)>LonLim(2)|| ...
        AEOLUS_latitude(j)<LatLim(1)||AEOLUS_latitude(j)>LatLim(2))
        AEOLUS_longitude(j) =[];
        AEOLUS_latitude(j)  =[];
        AEOLUS_altitude(j,:)=[];
        Extinction_355(j,:) =[];
    end
end

%筛选CALIPSO经纬度在研究区域的坐标点
Origin_Length=length(CALIPSO_longitude);
for i=1:Origin_Length
    j=Origin_Length-i+1;
    if(CALIPSO_longitude(j)<(LonLim(1)-0.3)||CALIPSO_longitude(j)>(LonLim(2)+0.3)|| ...
        CALIPSO_latitude(j)<(LatLim(1)-0.3)||CALIPSO_latitude(j)>(LatLim(2)+0.3))
        CALIPSO_longitude(j)=[];
        CALIPSO_latitude(j) =[];
        Backscatter_532(j,:)  =[];
        Extinction_532(j,:)   =[];
        Backscatter_1064(j,:) =[];
        Extinction_1064(j,:)  =[];
    end
end

clear i j Origin_Length;



%%
% 将CALIPSO的轨迹投影到AEOLUS的轨迹上
Fixed_CALIPSO_latitude=zeros([length(AEOLUS_latitude) , 1]);
Fixed_CALIPSO_longitude=zeros([length(AEOLUS_latitude) , 1]);
Fixed_Backscatter_1064=zeros([length(AEOLUS_latitude) , width(Backscatter_1064)]);
Fixed_Backscatter_532=zeros([length(AEOLUS_latitude) , width(Backscatter_1064)]);
Fixed_Extinction_1064=zeros([length(AEOLUS_latitude) , width(Backscatter_1064)]);
Fixed_Extinction_532=zeros([length(AEOLUS_latitude) , width(Backscatter_1064)]);


%记录该点离哪个点最近
Belong_flag=zeros(size(CALIPSO_latitude));
for i=1:length(CALIPSO_latitude)
    min_gap=1000;min_flag=0;
    for j=1:length(AEOLUS_latitude)
        temp_dif=abs(AEOLUS_latitude(j)-CALIPSO_latitude(i))...
        +abs(AEOLUS_longitude(j)-CALIPSO_longitude(i));
        if(temp_dif<min_gap)
            min_gap=temp_dif;
            min_flag=j;
        end
    end
    Belong_flag(i)=min_flag;
end
clear min_flag temp_dif min_gap ;

%按照AEOLUS的
Belong_count=zeros(size(AEOLUS_latitude));
for i=1:length(AEOLUS_latitude)
    %从属于AEOLUS第一个点
    temp_Lon=zeros(size(CALIPSO_longitude));
    temp_Lat=zeros(size(CALIPSO_latitude));
    temp_B532=zeros(size(Backscatter_532));
    temp_B1064=zeros(size(Backscatter_1064));
    temp_E532=zeros(size(Extinction_532));
    temp_E1064=zeros(size(Extinction_1064));
    for j=1:length(Belong_flag)
        if(Belong_flag(j)==i)%如果第j个CALIPSO确实从属第i个AEOLUS
            Belong_count(i)=Belong_count(i)+1;%AEOLUS的从属点加1
            temp_Lon(j)=CALIPSO_longitude(j);
            temp_Lat(j)=CALIPSO_latitude(j);
            temp_B532(j,:)=Backscatter_532(j,:);
            temp_B1064(j,:)=Backscatter_1064(j,:);
            temp_E532(j,:)=Extinction_532(j,:);
            temp_E1064(j,:)=Extinction_1064(j,:);
        end
    end
    %根据temp的内容计算平均值
    if(Belong_count(i)==0)
        continue;
    else
        %阅历temp根据有效值数量分别计算各参数
        temp_count=0;temp_total_lon=0;temp_total_lat=0;
        for j=1:length(temp_Lon)
            if(temp_Lat(j)~=0&&temp_Lon(j)~=0)
                temp_count=temp_count+1;
                temp_total_lon=temp_total_lon+temp_Lon(j);
                temp_total_lat=temp_total_lat+temp_Lat(j);
            end
        end
        %检查属于AEOLUS(i)的点是否计算有误，无误则计算每个剖面的Backscatter和Extinction
        if(Belong_count(i)~=temp_count)
            disp('wrong with the temp');
            pause;
        else
            %1完成经纬度的resize
            Fixed_CALIPSO_latitude(i)=temp_total_lat/temp_count;
            Fixed_CALIPSO_longitude(i)=temp_total_lon/temp_count;
            clear temp_total_lon temp_total_lat;
            
            %开始计算从属于第i个AEOLUS的 Backscatter和Extinction
            %2B532
            for j=1:width(temp_B532) %j代表高度层,399个
               temp_count=0;temp_total=0;
               for k=1: length(temp_B532)%k代表行,即哪个CALIPSO点
                   if(temp_B532(k,j)~=0.0&&temp_B532(k,j)~=-9999)
                       temp_count=temp_count+1;
                       temp_total=temp_total+temp_B532(k,j);
                   end
               end
               %temp_count=0时代表该层属于AEOLUS(I)无有效数据,将被置为NaN
                Fixed_Backscatter_532(i,j)=temp_total/temp_count;
            end

            %3B1064
            for j=1:width(temp_B1064) %j代表高度层,399个
               temp_count=0;temp_total=0;
               for k=1: length(temp_B1064)%k代表行,即哪个CALIPSO点
                   if(temp_B1064(k,j)~=0.0&&temp_B1064(k,j)~=-9999)
                       temp_count=temp_count+1;
                       temp_total=temp_total+temp_B1064(k,j);
                   end
               end
               %temp_count=0时代表该层属于AEOLUS(I)无有效数据,将被置为NaN
                Fixed_Backscatter_1064(i,j)=temp_total/temp_count;
            end
            
            %4 E532
            for j=1:width(temp_E532) %j代表高度层,399个
               temp_count=0;temp_total=0;
               for k=1: length(temp_E532)%k代表行,即哪个CALIPSO点
                   if(temp_E532(k,j)~=0.0&&temp_E532(k,j)~=-9999)
                       temp_count=temp_count+1;
                       temp_total=temp_total+temp_E532(k,j);
                   end
               end
               %temp_count=0时代表该层属于AEOLUS(I)无有效数据,将被置为NaN
                Fixed_Extinction_532(i,j)=temp_total/temp_count;
            end
           
            %5 E1064
            for j=1:width(temp_E1064) %j代表高度层,399个
               temp_count=0;temp_total=0;
               for k=1: length(temp_E1064)%k代表行,即哪个CALIPSO点
                   if(temp_E1064(k,j)~=0.0&&temp_E1064(k,j)~=-9999)
                       temp_count=temp_count+1;
                       temp_total=temp_total+temp_E1064(k,j);
                   end
               end
               %temp_count=0时代表该层属于AEOLUS(I)无有效数据,将被置为NaN
                Fixed_Extinction_1064(i,j)=temp_total/temp_count;
            end

        end
    end

end


%删除没有匹配点的AEOLUS同步删除空CALIPSO
Origin_Length=length(Belong_count);
for i=1:Origin_Length
    j=Origin_Length-i+1;
    if(Belong_count(j)==0)%没有匹配的CALIPSO
    Fixed_CALIPSO_latitude(j)  =[];
    Fixed_CALIPSO_longitude(j) =[];
    AEOLUS_longitude(j)        =[];
    AEOLUS_latitude(j)         =[];
    AEOLUS_altitude(j,:)=[];
    Fixed_Backscatter_1064(j,:)=[];
    Fixed_Backscatter_532(j,:) =[];
    Fixed_Extinction_1064(j,:) =[];
    Fixed_Extinction_532(j,:)  =[];
    Extinction_355(j,:)        =[];
    end
end
clear temp_B532 temp_B1064 temp_E532 temp_E1064 temp_Lat temp_Lon ;
clear temp_total temp_count Origin_Length i j k  Extinction_1064 Extinction_532;
clear  Backscatter_1064 Backscatter_532 CALIPSO_latitude CALIPSO_longitude;
clear Belong_count Belong_flag;


%无效值记为NaN

Extinction_355(Extinction_355==-1000000)=NaN;
Extinction_355(Extinction_355==0)=NaN;



%最后结果
CALIPSO_latitude=Fixed_CALIPSO_latitude;
CALIPSO_longitude=Fixed_CALIPSO_longitude;
CALIPSO_altitude=zeros(1,24);
Backscatter_1064=zeros(length(AEOLUS_latitude),24);
Backscatter_532=zeros(length(AEOLUS_latitude),24);
Extinction_1064=zeros(length(AEOLUS_latitude),24);
Extinction_532=zeros(length(AEOLUS_latitude),24);

%CALIPSO的高度层同化 
%7*2 8*5 21*5 20*12 
%i代表每个合成后的footprint
%j代表每个高度层


for j=1:24  %7*2
    %1 altitude 均有效
    Floor_high=15;
    temp_count=0;temp_total=0;
    %k具体参与合成的原高度层
    for k=(j-1)*Floor_high+40 :(j-1)*Floor_high+Floor_high+39
        temp_count=temp_count+1;
        temp_total=temp_total+Fixed_CALIPSO_altitude(k);
    end
    CALIPSO_altitude(j)=temp_total/temp_count;

    for i=1:length(AEOLUS_latitude)
        %2 Backscatter_532
        temp_count=0;temp_total=0;
        for k=(j-1)*Floor_high+1 :(j-1)*Floor_high+Floor_high
            if(~isnan(Fixed_Backscatter_532(i,k)))
                temp_count=temp_count+1;
                temp_total=temp_total+Fixed_Backscatter_532(i,k);
            end
        end
        Backscatter_532(i,j)=temp_total/temp_count;
        %3 Backscatter_1064
        temp_count=0;temp_total=0;
        for k=(j-1)*Floor_high+1 :(j-1)*Floor_high+Floor_high
            if(~isnan(Fixed_Backscatter_1064(i,k)))
                temp_count=temp_count+1;
                temp_total=temp_total+Fixed_Backscatter_1064(i,k);
            end
        end
        Backscatter_1064(i,j)=temp_total/temp_count;
        %4 Extinction_532
        temp_count=0;temp_total=0;
        for k=(j-1)*Floor_high+1 :(j-1)*Floor_high+Floor_high
            if(~isnan(Fixed_Extinction_532(i,k)))
                temp_count=temp_count+1;
                temp_total=temp_total+Fixed_Extinction_532(i,k);
            end
        end
        Extinction_532(i,j)=temp_total/temp_count;
        %5 Extinction_1064
        temp_count=0;temp_total=0;
        for k=(j-1)*Floor_high+1 :(j-1)*Floor_high+Floor_high
            if(~isnan(Fixed_Extinction_1064(i,k)))
                temp_count=temp_count+1;
                temp_total=temp_total+Fixed_Extinction_1064(i,k);
            end
        end
        Extinction_1064(i,j)=temp_total/temp_count;
    end

 
end




%
clear Fixed_Extinction_1064 Fixed_Extinction_532 Fixed_Backscatter_1064 Fixed_Backscatter_532;
clear Fixed_CALIPSO_latitude Fixed_CALIPSO_longitude i j k Floor_high temp_total temp_count;


end

2.使用案例

代码如下（示例）：

file_name_AEOLUS='AE_OPER_ALD_U_N_2A_20200914T225411022_005412007_011957_0004.DBL';
file_name_CALIPSO='CAL_LID_L2_05kmAPro-Standard-V4-21.2020-09-14T18-39-46ZD.hdf';
file_path_CALIPSO=[root_path_CALIPSO,file_name_CALIPSO];
[AEOLUS_altitude,AEOLUS_latitude,AEOLUS_longitude,...
   CALIPSO_altitude, CALIPSO_latitude,CALIPSO_longitude,...
   Backscatter_532,Backscatter_1064,Extinction_355,Extinction_532,...
   Extinction_1064] = resize(file_name_CALIPSO,file_name_AEOLUS);

---