CALIPSO数据下载与MATLAB读取

目录

官网请求数据

下载方法

MATLAB导入数据

level 2 产品说明+代码：以Feature_Classification_Flags为例

matlab代码重现

level 1 matlab代码

---

官网请求数据

登陆官网，选择自己需要的数据产品、时间、范围等

conform request

几分钟后（有一次是几个小时后）邮箱中会有下载邮件。另外需要下载、安装wget

据说wget的安装后还需要添加环境变量，不然可能会闪退。（本人没加....）

可以直接在设置中查找"高级系统设置”

 双击path加入安装地址。

下载方法

打开cmd，输入要进入的磁盘，如进入E盘：C:\Users\Administrator>e: 【注意加冒号】

然后进入（cd）具体存放数据的文件夹：

输入下图邮件中的字符+链接

MATLAB导入数据

FILE_NAME_L1 = ['E:\PBLH\Code\originalData\CALIPSO\', ...
'2022070004946_65124\CAL_LID_L1-Standard-V4-10.2019-02-14T04-29-30ZD_Subset.hdf']; %文件地址

info_L1 = hdfinfo(FILE_NAME_L1);%读取hdf文件信息
dsetsL1 = info_L1.SDS;%读取SDS信息，包括具体变量的名字
% 读取具体变量
UTC_Time = hdfread(FILE_NAME_L1, 'Profile_UTC_Time');% 读取UCT时间
Time= hdfread(FILE_NAME_L1, 'Profile_Time');% TAI时间，用秒记录 
TAB = hdfread(FILE_NAME_L1, 'Total_Attenuated_Backscatter_532');%
PAB = hdfread(FILE_NAME_L1, 'Perpendicular_Attenuated_Backscatter_532');%
lat = hdfread(FILE_NAME_L1, 'Latitude');%纬度
lon = hdfread(FILE_NAME_L1, 'Longitude');%经度
……

level 2 产品说明+代码：以Feature_Classification_Flags为例

数据结构如图所示，代表分类结果的变量 Feature_Classification_Flags 维度为 (3712, 5515)，3712 是沿卫星轨道扫过的点数，5515 是每个点的垂直结构——本来是一个二维的数组，现展平成长为 5515 的一维数组，在使用时需要恢复它的二维排布。

官网（链接如下）说明，二维数组分为垂直三层，每一层的水平分辨率和垂直分辨率都不一样。二维数组沿轨道方向长 5 km，高 30 km，文件中自带的经纬度和时间对应 5 km 中点处的值。实际使用时为了方便起见，取每一层的第一条廓线来代表此层的结果，也就是图中标红的那三条。于是每个点的二维数组简化为垂直分辨率不一的一条廓线。CALIPSO - User’s Guide

FROM: Python 处理CALIPSO L2 VFM激光雷达数据及绘图【以3.15北方沙尘暴为例】

来自官网的综合实例：Comprehensive Examples of LaRC Products

以及官网matlab level2关于'Feature_Classification_Flags'的实例：Level 2 Vertical Feature Mask Version 3 file in MATLAB. 

3 profiles for 20.2km (base) to 30.1km (top) @ 180m

 (index 1-165 / 55 samples per profile)

即最上层层是3条廓线：

3 profiles mean horizontal resolution is 1667m because 3 * 1667m = 5km.

55 samples mean vertical resolution is 180m because 55 * 180m = 9.9km  = 30.1km - 20.2km.

5 profiles for 8.2km (base) to 20.2km (top) @ 60m

 (index 166 - 1165 / 200 samples per profile)

15 profiles for -0.5km (base) to 8.2km (top) @ 30m 

(index 1166 - 5515 / 290 samples per profile)

matlab代码重现

最主要的是高度的推算，注意数据的转置。

%% L2
FILE_NAME_L2 = ['E:\PBLH\Code\originalData\CALIPSO\', ...
'2022061143028_64875\CAL_LID_L2_VFM-ValStage1-V3-41.2021-11-29T23-32-39ZD.hdf'];
info_L2 = hdfinfo(FILE_NAME_L2);%读取hdf文件信息
dsets = info_L2.SDS;%读取SDS信息
lat = hdfread(FILE_NAME_L2, 'Latitude');%纬度
lon = hdfread(FILE_NAME_L2, 'Longitude');%经度
data = hdfread(FILE_NAME_L2, 'Feature_Classification_Flags');%

% Select the 1-3 bits for Feature Type data. 
data = bitand(data, 7);
% Subset latitude values for the region of interest (40N to 62N).
lat = lat(3501:4000);
dim = size(lat, 1);
% -0.5km to 8.2km
data2d = squeeze(data(3501:4000, 1166:5515));   
data3d = reshape(data2d,[dim, 290, 15]); %(8.2+0.5)/0.03=290
data = squeeze(data3d(:,:,1));
% Focus on cloud (=2) data only.
data(data > 2) = 1;
data(data < 2) = 1;
% -0.5km to 8.2km
for i=0:289
 altitude(i+1) = -0.5 + i*0.03;
end 
% Create a Figure to Plot the data.
f = figure('Name', FILE_NAME_L2, ...
    'Renderer', 'zbuffer', ...
    'Position', [0,0,800,600], ...
    'Visible', 'off', ...
    'PaperPositionMode', 'auto');
% Create a custom color map for 2 different Feature Type key values.
% 为两个不同的要素类型键值创建自定义颜色贴图。
cmap=[                       %  Key            R   G   B
      [0.00 0.00 1.00];  ... %  1=not cloud   [000,000,255]
      [1.00 1.00 1.00];  ... %  2=cloud       [255,255,255]
     ];     

colormap(cmap);
caxis([1 2]); 
contourf(lat, altitude, rot90(data, 1));
% Set axis labels.
xlabel('Latitude (degrees north)'); 
ylabel('Altitude (km)');
% Put colorbar.
y = [1, 2];
h = colorbar('YTickLabel', {'Others', 'Cloud'}, 'YTick', y);
tstring = {FILE_NAME_L2;['Feature Type  (Bits 1-3) in' ...
                    ' Feature Classification Flag']};
title(tstring, 'Interpreter', 'none', 'FontSize', 16, ...
      'FontWeight','bold');

hold off;
saveas(f, [FILE_NAME_L2 '.v.m.png']);
exit;

level 1 matlab代码

'Total_Attenuated_Backscatter_532'

FILE_NAME_L1 = ['E:\PBLH\Code\originalData\CALIPSO\', ...
'2022070004946_65124\CAL_LID_L1-Standard-V4-10.2019-02-14T04-29-30ZD_Subset.hdf'];
info_L1 = hdfinfo(FILE_NAME_L1);%读取hdf文件信息
dsetsL1 = info_L1.SDS;%读取SDS信息
UTC_Time = hdfread(FILE_NAME_L1, 'Profile_UTC_Time');%
TAB = hdfread(FILE_NAME_L1, 'Total_Attenuated_Backscatter_532');%
PAB = hdfread(FILE_NAME_L1, 'Perpendicular_Attenuated_Backscatter_532');%
lat = hdfread(FILE_NAME_L1, 'Latitude');%纬度
lon = hdfread(FILE_NAME_L1, 'Longitude');%经度

% -0.5km to 8.3km
TAB2d = TAB(:, 289:578)';   
PAB2d = PAB(:, 289:578)';   
for i=1:size(lon,1) 
    TAB3d(:,i) = SLidar_SG_filter(TAB2d(:,i)); 
end
for i=0:289
 altitude(i+1) = -0.5 + i*0.03;%(8.3+0.5)/0.03=293.3
end 

figure
set(gca,'Fontsize',12,'FontWeight','bold','Fontname','Times New Ramon','color','none','linewidth',1.5);
h=mesh(lat,altitude,TAB3d);
axis xy; axis tight;
view(0,90);colormap(jet);
colorbar
ylabel('Altitude [km]','Fontsize',12,'FontWeight','bold','Fontname','Times New Roman');
xlabel('Latitude','Fontsize',12,'FontWeight','bold','Fontname','Times New Roman');