depolarization_ratio_frequency.rar_calipso_气溶胶_退偏比

clear all %清除工作区所有变量 clc %清除命令窗口命令记录 path_data='G:\calipso\level2\alay\2016\1\'; % 将字符串（数据所在文件夹路径）赋值给path_data files0=dir(path_data); % dir函数，得到path_data路径下，所有文件和文件夹的名称，赋值给变量files0 files=files0(arrayfun(@(x) ~strcmp(x.name(1),'.'),files0)); % 将首字符为“.”的文件或文件夹剔除，注：首字符为“.”的为隐藏文件（夹） files=files(arrayfun(@(x) strcmp(x.name(1:3),'CAL'),files)); % 获得前三个字符为CAL的所有文件名称 nfiles=length(files); % 文件个数 i=0; for a=0:5 for b=0:1 for c=0:3 for d=0:3 for e=0:3 i=i+1; allvalue(i)=a*8192+b*4096+c*128+d*32+e*8; % Feature_Classification_Flags 中可以判断为某一类气溶胶（需加上特定数值）所有可能数值 end end end end end count1=0; % 计数 count2=0; count3=0; count4=0; count5=0; count6=0; count7=0; for k=1:nfiles % 循环，每个数据文件都统计 dirname=files(k).name; %第k个数据文件的文件名 Latitude = hdfread([path_data,dirname], '/Latitude'); %读取数据文件中的纬度变量 Longitude = hdfread([path_data,dirname], '/Longitude'); % 读取数据文件中的经度变量 Feature_Classification_Flags = hdfread([path_data,dirname], '/Feature_Classification_Flags'); % 读取 Feature Classification Flags 数据，可根据此数据判断气溶胶类型 Integrated_Particulate_Depolarization_Ratio = hdfread([path_data,dirname], '/Integrated_Particulate_Depolarization_Ratio'); % 读取层积分退偏比数据 % 类似，若读取色比数据，则用以下命令： % Integrated_Particulate_Color_Ratio = hdfread([path_data,dirname], '/Integrated_Particulate_Color_Ratio'); % 注：退偏比取值范围一般为0-1.0之间，色比取值范围要更大，若要统计色比数据，需合稍加修改程序 Lat=Latitude(:,2); % 中心纬度 Lon=Longitude(:,2); % 中心经度 flagx=find(Lon>70 & Lon <105); % 找出经度 >70° 且 <105°的数据 % 根据实际需要修改 flagy=find(Lat>27 & Lat <40); % 找出纬度 >27° 且 <40°的数据 % 根据实际需要修改 domain=intersect(flagx,flagy); % intersect函数，找出两个数组的交集，这里即找出 70°<Lon<105° & 27°<Lat<40°的数据 FCF=Feature_Classification_Flags(domain,:); % 得到满足经纬度范围条件的 Feature Classification Flags 数据 PDR=Integrated_Particulate_Depolarization_Ratio(domain,:); % 得到满足经纬度范围条件的退偏比数据 flag=find(ismember(FCF,allvalue+1027)==1); % 找出气溶胶类型为沙尘的数据记录 % allvalue+1027 即为沙尘气溶胶的所有可能的Feature_Classification_Flags值 % allvalue+515 则气溶胶类型为 clean marine % allvalue+1539 则气溶胶类型为 polluted continental % allvalue+2051 则气溶胶类型为 clean continental % allvalue+2563 则气溶胶类型为 polluted dust % allvalue+3075 则气溶胶类型为 smoke % 因此可以用此命令统计不同类型气溶胶的频率分布 % flag=find(ismember(FCF,allvalue+ X )==1)得到某气溶胶的数据记录 % count=length(flag); 即可得到记录次数 % ismember函数，用法： C=ismember(A,B); 判断A数组中的数是否等于B数组中的任意一个数，返回一个和A等大小的逻辑数组C % 例如：A=[1,2,3;4,5,6]; B=[2,8;6,9] % C=ismember(A,B); 则 C=[0,1,0;0,0,1]; class1=find(PDR>0 & PDR<=0.1); % 找出退偏比 >0 且 ≤0.1 的退偏比的数据记录 class2=find(PDR>0.1 & PDR<=0.2); % 找出退偏比 >0.1 且 ≤0.2 的退偏比的数据记录 class3=find(PDR>0.2 & PDR<=0.3); % 找出退偏比 >0.2 且 ≤0.3 的退偏比的数据记录 class4=find(PDR>0.3 & PDR<=0.4); % 找出退偏比 >0.3 且 ≤0.4 的退偏比的数据记录 class5=find(PDR>0.4 & PDR<=0.5); % 找出退偏比 >0.4 且 ≤0.5 的退偏比的数据记录 class6=find(PDR>0.5 & PDR<=0.6); % 找出退偏比 >0.5 且 ≤0.6 的退偏比的数据记录 class7=find(PDR>0.6 & PDR <=1.0); % 找出退偏比 >0.6 且 ≤1.0 的退偏比的数据记录 count1=count1+length(intersect(flag,class1)); % 求和，退偏比范围在0-0.1的数据记录次数，下面类似 count2=count2+length(intersect(flag,class2)); count3=count3+length(intersect(flag,class3)); count4=count4+length(intersect(flag,class4)); count5=count5+length(intersect(flag,class5)); count6=count6+length(intersect(flag,class6)); count7=count7+length(intersect(flag,class7)); end total=count1+count2+count3+count4+count5+count6+count7; % 数据记录总数 f1=count1/total; % 沙尘气溶胶退偏比范围为0-0.1的出现频率 f2=count2/total; f3=count3/total; f4=count4/total; f5=count5/total; f6=count6/total; f7=count7/total;