【calipso廓线数据下载&处理】

**

calipso气溶胶廓线数据长时间统计分析

**

calipso数据下载

calispo数据下载地址已经更新，网址如下：https://subset.larc.nasa.gov/calipso/login.php
需要注册NASA数据下载账号，建议使用gmail邮箱注册（科学上网）。
1）.选择数据内容
2）.选择时间范围
3）.选择空间范围
4）.确认订单，提交订单
完成以上步骤后，数据下载链接会发送到你的邮箱中，下载链接是ftp下载方式，建议使用火狐/Google游览器打开，并在浏览器内安装一下down them all插件！！！按照文件类型筛选下载链接，再一键下载。

以消光、退偏、后向散射系数为例，统计分析2015-2020年的三维空间分布特征为例
数据：L2 apro数据

主要流程：

1.数据筛选，跟自己的需求筛选数据
筛选出指定年份/季节/月份的数据，考虑到信噪比问题，尽量选择夜间数据（*ZN_Subset.hdf）
2.以指定分辨率建立格网，格网size=nl x ns x 399
nl–>行数
ns–>列数
399–>层数，calipso L2 apro数据为399层
3.对筛选出的文件名进行for循环，按照经纬度信息，将光学参数填入格网（累计），并记录每个格网的填入的点数目。
4.计算光学参数的均值
光学参数均值=累计值/像元数目

注意：
在步骤3中，需要根据CAD（质量标识码进行数据质量控制，CAD<0为气溶胶，反之则为非气溶胶，CAD的绝对值代表置信度，一般选择CAD<-20作为气溶胶质控标准）；VAD作为分类标识码，在研究指定类型气溶胶的时候，再结合VAD信息，进行指定类型气溶胶的筛选。

下面开始数据处理：

%% 只统计了气溶胶，未细分具体类型
clear all;
close all;
clc;
%% 以全部数据为例（选择较低信噪比的夜间观测资料）
rootdir = '/desktop/L2Apro/';
outpath = '/desktop/L2Apro/out/';
[filename_list] = get_filen(rootdir,'ZN_Subset.hdf');
%% 创建空数组
lonst = 0 
loned = 180
latst = 0
lated = 90
res = 0.5
nl = int((lated-latst)/res);
ns = int((loned-lonst)/res);
ec_sum = zeros(nl,ns,399);
tbc_sum = zeros(nl,ns,399);
pdr_sum = zeros(nl,ns,399);
pixel_sum = = zeros(nl,ns,399);
for i = 1:size(filename_list,1)
	filename = strcat(rootdir,filename_list(i,:));
	[tbc,pdr,ec,lon,lat] = read_apro(filename);
	num = length(lon);
	for j = 1:num
		xx = round((lon(j)-lonst)/res)+1;
		yy = round((lat(j)-latst)/res)+1;
		if xx>0 && xx<=nl && yy>0 && yy<=ns 
			for h=1:399
				ec_sum(xx,yy,h) = ec_sum(xx,yy,h) + ec(j,h);
				pdr_sum(xx,yy,h) = pdr_sum(xx,yy,h) + ec(j,h);
				tbc_sum(xx,yy,h) = tbc_sum(xx,yy,h) + ec(j,h);
				pixel_sum(xx,yy,h) = pixel_sum(xx,yy,h) + 1;
			end
		end
	end	
end
%% 计算光学参数的均值
pixel_sum(pixel_sum==0) = 1;
ec_m = ec_sum./pixel_sum;
pdr_m = pdr_sum./pixel_sum;
tbc_m = tbc_sum./pixel_sum;
ec_m(ec_m==0) = nan;
pdr_m(pdr_m==0) = nan;
tbc_m(tbc_m==0) = nan;
save(strcat(outpath,'ec.mat'),'ec_m');
save(strcat(outpath,'pdr.mat'),'pdr_m');
save(strcat(outpath,'tbc.mat'),'tbc_m');



%%
function [tbc,pdr,ec,lon,lat] = read_apro(filename)
	metadata = hdfread(strname, '/metadata', 'Fields', 					'Lidar_Data_Altitudes', 'FirstRecord',1 ,'NumRecords',1);
    alts = metadata{1};
    %%
    variable = 'Latitude';
    lat = hdfread(strname,variable);
	variable = 'Longitude';
    lon = hdfread(strname,variable);
    attribute = 'Total_Backscatter_Coefficient_532';
    TBC = hdfread(strname,attribute);
    TBC = TBC';

    attribute = 'Particulate_Depolarization_Ratio_Profile_532';
    PDR = hdfread(strname,attribute);
    PDR = PDR';

    attribute = 'Surface_Elevation_Statistics';
    SES = hdfread(strname,attribute);
    SES = SES(:,3)';

    attribute = 'Backscatter_Coefficient_1064';
    BC = hdfread(strname,attribute);   
    BC = BC';
                    %---
    attribute = 'Extinction_Coefficient_532';
    EC = hdfread(strname,attribute);
    EC = EC';

    attribute = 'Extinction_Coefficient_Uncertainty_532';
    ECU = hdfread(strname,attribute);
    ECU = ECU';

    attribute = 'CAD_Score';
    CAD = hdfread(strname,attribute);
    CAD3 = CAD(:,1:399,:);
    CAD1 = (CAD3(:,:,1:1))';
    CAD2 = (CAD3(:,:,2:2))';

    attribute = 'Atmospheric_Volume_Description';
    AVD = hdfread(strname,attribute);
    AVD3 = AVD(:,1:399,:);
    AVD1 = (AVD3(:,:,1:1))';
    AVD2 = (AVD3(:,:,2:2))';
    %%
    ftype1 = bitand(AVD1, 7);
    ftype2 = bitand(AVD2, 7);
    typeBits1 = bitshift(AVD1, -9);
    typeBits2 = bitshift(AVD2, -9);
    subtype1 = bitand(typeBits1, 7);
    subtype2 = bitand(typeBits2, 7);
    %%
    [nn,mm] = size(BC);
    useSamples = true(nn,mm);
    gg_d = ones(nn,mm);
    useSamples(EC == -9999) = false;

    [ftype, subtype] = get_feature_type_and_subtype(AVD);

    isAerosol = squeeze(ftype(1,:,:) == 3 | ftype(2,:,:) == 3);
    useSamples(~isAerosol) = false;

    goodCAD = squeeze( CAD1 < -20 & ...
                       CAD2 < -20 );
    gg_d = gg_d.*goodCAD;

    useSamples(~goodCAD) = false;
	indx = useSamples.*gg_d;
	ec = EC.*indx;
	pdr = PDR.*indx;
	tbc = TBC.*indx;

end

function [filename_list] = get_filen(inpath,attrs):
	%% inpath:输入路径
	%% attrs:数据类型
	File = dir(fullfile(inputdir,'*.hdf'));%遍历文件，D是白天，N	是夜晚，数据是APro数据
	strtxt = {File.name}';
	filename_list = char(strtxt);
end


function [ftype, subtype] = get_feature_type_and_subtype(FCF)
%此函数接受功能分类标志(FCF)数组，返回数组功能类型和特征子类型。
%功能和子类型在CALIPSO数据产品目录(版本3.2)的表44中定义。
%请注意，要素子类型的定义取决于要素类型
%输入FCF，一组CALIPSO功能分类标志。Type：uint16%输出%ftype，与包含该功能类型
%的FCF大小相同的数组。Type：uint16子类型，与包含特征子类型的FCF大小相同的数组。

% ---------- 获取特征类型 --------- %
% 这些是返回值及其定义
%   0 = invalid (bad or missing data)
%   1 = "clear air"
%   2 = cloud
%   3 = aerosol
%   4 = stratospheric feature; polar stratospheric cloud (PSC) or stratospheric aerosol %平流层特性;极地平流层云(PSC)或平流层气溶胶
%   5 = surface地表
%   6 = subsurface 地表下
%   7 = no signal (totally attenuated)%无信号(完全衰减)
    ftype = bitand(FCF, 7);
    

%     ftype = bitand(FCF, 7);

% -------- Get feature subtype -------- % 
% 备注：特征子类型的定义将取决于功能类型！有关详细信息，请参阅CALIPSO Data%Products Catalog中的功能分类标志定义表。

% -- 气溶胶类型定义
%   0 = not determined
%   1 = clean marine
%   2 = dust
%   3 = polluted continental
%   4 = clean continental
%   5 = polluted dust
%   6 = smoke
%   7 = other

% -- 云类型定义
%   0 = low overcast, transparent
%   1 = low overcast, opaque
%   2 = transition stratocumulus
%   3 = low, broken cumulus
%   4 = altocumulus (transparent)
%   5 = altostratus (opaque)
%   6 = cirrus (transparent)
%   7 = deep convective (opaque)%深对流(不透明)

   typeBits       = bitshift(FCF, -9);
    subtype = bitand(typeBits, 7);

end

后续再更新绘图代码

数据处理后，绘图