提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档
前言
本文用于记录分享处理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);