HDF4或者HDF-EOS读取与可视化

最新推荐文章于 2024-07-06 18:37:49 发布

大作家佚名

最新推荐文章于 2024-07-06 18:37:49 发布

阅读量8.4k

点赞数 11

分类专栏： MATLAB 算法/方法 CALIPSO

本文链接：https://blog.csdn.net/wokaowokaowokao12345/article/details/79115443

版权

MATLAB 同时被 3 个专栏收录

48 篇文章 54 订阅

订阅专栏

算法/方法

37 篇文章 6 订阅

订阅专栏

CALIPSO

16 篇文章 44 订阅

订阅专栏

HDF4或者HDF-EOS读取

引言

这里介绍如何使用Maltlab自带函数读取HDF4文件中的数据，实例数据文件是CALIPSO获取的气溶胶廓线数据。

HDF4

CALIPSO，CloudSat和MODIS选择HDF4（Hierarchical Data Format）作为存储和分发科学数据和辅助数据的主要格式。它是一种自我描述的，与平台无关的二进制格式，用于存储大型数据集和元数据。可以包含在HDF文件中的对象包括多维数组（Scientic Data Set 或SDS）、表数据，即元数据（Vdata）、栅格图像、调色板和注释。每个对象以及本身可以进一步用任意数量的名称 -值对（’属性’）来描述。

HDF数据集类型

当前的HDF数据集类型包括3种类型：
1. Scientific data sets (SDS)科学数据集，包含多维矩阵，和数据的描述信息矩阵；
2. Raster image set(RIS)栅格数据集，包括栅格图像和相应的描述信息；
3. Vset(VS)，包含用户希望包含的任何种类的HDF对象的通用分组结构。

函数语法

这里介绍的hdf读取函数针对HDF 4.x版本，或者HDF-EOS 2.x版本的文件。如果要读取HDF5，则需要使用h5read函数。

data = hdfread(filename, datasetname)
data = hdfread(hinfo)
data = hdfread(...,param,value,...)
data = hdfread(filename,EOSname,param,value,...)
[data,map] = hdfread(...)

函数参数讲解

data = hdfread（filename，datasetname）
返回由filename指定的HDF4或HDF-EOS文件中datasetname指定的数据集中的所有数据。可以使用hdfinfo确定HDF4中的数据集名称。
data = hdfread（hinfo）返回由hdfinfo函数返回的structurehinfo指定的数据集中的所有数据。指定hinfo结构中与特定类型数据集相关的字段，并使用索引来指定哪个数据集有多个数据集。请参阅指定要读取的数据集以获取更多信息。
data = hdfread(…,param,value,…)根据指定的参数和值对返回数据的子集。
data = hdfread(filename,EOSname,param,value,…)返回由EOSname，参数和值指定的HDF-EOS点，格网或条（与卫星轨道垂直扫描的条带）数据字段。
[data,map] = hdfread(…)返回栅格数据的影像图和颜色条。

CALIPSO数据产品

这里只关注CALIPSO中的CALIOP激光雷达的数据产品，下图是CALIPSO数据产品中的命名规则。

CALIPSO数据命名规则

这篇博客里使用的是气溶胶数据，数据为“CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf”，根据上图命名规则，我们可以知道，这个数据是CALIPSO激光雷达CALIOP的Level2 5km分辨率的气溶胶廓线产品，版本是V3，发布时间是2010年4月17日20时9分13秒白天的观测数据。知道了这个我们就可以更有针对性的对HDF数据进行读取理解了。

关于矩阵大小

消光系数的二维矩阵及其大小为：numberAlts×numberProfiles，分别表示高度和廓线数量。CAD Score，Extinction QC和大气体积描述等描述性标志都是三维矩阵，大小为：2×numberAlts×numberProfiles。

实例

代码

clc
clear all
fileinfo = hdfinfo('CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf')

结果

fileinfo = 

      Filename: '...CAL_LID_L2_05kmAPro-Prov-V3-01.2...'
    Attributes: [1x2 struct]
           SDS: [1x47 struct]
         Vdata: [1x1 struct]

科学数据集（SDS）提供了一个框架，用于存储多维数据阵列，并提供增强数据的描述性信息。当前规格支持SDS阵列中的以下类型的数字。
下表总结了激光雷达廓线产品中包含的每个科学数据集（SDS）的内容。每个参数都使用相应HDF文件中使用的相同SDS名称列出。

NASA官方文档说明
这里写图片描述

Matlab实际读取后的SDS信息

激光雷达气溶胶廓线数据产品包括三个Vdata记录类型（即元数据），激光雷达2级气溶胶廓线数据特有的元数据参数产品列于下表。

NASA官方文档说明
这里写图片描述

Matlab实际读取后的Vdata信息

读取SDS数据集中的温度信息

sds_info = fileinfo.SDS(21)

结果

sds_info = 

       Filename: '...CAL_LID_L2_05kmAPro-Prov-V3-0...'
           Type: 'Scientific Data Set'
           Name: 'Temperature'
           Rank: 2
       DataType: 'float'
     Attributes: [1x3 struct]
           Dims: [2x1 struct]
          Label: {}
    Description: {}
          Index: 20

读取温度的三种方法
方法1

Temperature = hdfread(sds_info);

方法2
指定文件，字段，默认输出全部信息

Temperature = hdfread('CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf', 'Temperature');

方法3
指定文件，字段，索引

Temperature = hdfread('CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf', 'Temperature', 'Index', {[1  1],[1  1],[3744   399]});

该数据是3744X399的二维矩阵，‘Index’的变量意义是：｛开始，步长，读取个数｝，由于是二维矩阵，所以分别指定了两维的起始、步长和个数。

提取CALIPSO数据产品所需Matlab函数总结

permute

该函数具有矩阵转置、变形的功能，比如A矩阵是一个NXMXD的三维矩阵，
语法：B = permute(A,[3 2 1]); 的作用就是将A转换成DXMXN的B矩阵。
[3 2 1]就是将原先的第3维变换为新矩阵的第一维，第二维还是第二维，第一维变成第3维。

dec2bin 和bin2dec

十进制和二进制互转函数，语法：a = dec2bin(b)

bitand

返回两个数值型数值在按位进行AND运算后的结果。
语法：C = bitand(A,B)
实例
x = 5 && 二进制为 0101
y = 6 && 二进制为 0110
bitand(x,y) && 返回值 4，二进制为 0100

bitand(5,6)

结果4

bitshift

返回向左移动给定位后的结果。语法：C = bitshift(A,K) 。
实例
x = 5 &&　二进制为 0101
y = 1 &&　左移一位
bitlshift(x,y) && 返回值 10，二进制为 1010

bitshift(5,1)

结果10

squeeze

移除矩阵单一维。
语法：B = squeeze(A)
实例
随机产生一个1x2x3的矩阵A，然后squeeze（A）将返回一个2x3的矩阵，将第一维去掉(因为第一维大小为1)。

A = rand(1,2,3)
B = squeeze(A)

A(:,:,1) = 0.8147    0.9058
A(:,:,2) = 0.1270    0.9134
A(:,:,3) = 0.6324    0.0975
B =
    0.8147    0.1270    0.6324
    0.9058    0.9134    0.0975

逻辑数组（矩阵）

Matlab逻辑类型的数据只有逻辑真和逻辑假，用1表示逻辑真，用函数true()表示；用0表示逻辑假，用函数false()表示。在Matlab中，可以用true()和false()函数创建逻辑矩阵。
语法

true(m)
true(m,n)

应用于从矩阵（数组）筛选感兴趣数据，实例

A = [2 3 54 123 48 13 5 7];
a = [1 0 1 0 0 0 1 0];
b = logical(a)；
A(b)

结果2 54 5

clc
clear all
close all

% 读取文件信息
fileinfo = hdfinfo('CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf');

% 读取高度数据
metadata = hdfread(fileinfo.Filename, '/metadata', 'Fields', 'Lidar_Data_Altitudes', 'FirstRecord',1 ,'NumRecords',1);
alts = metadata{1};
numberAlts = length(alts);

%% 基于readHDF读取矩阵信息，第一个元素索引是0，不同于Matlab矩阵索引
% 这个例子中，我们指定纬度范围，并提取该范围的廓线信息
latMin = 15;
latMax = 25;

%% 读取纬度信息
product_name = 'Latitude';
[info, lat] = readHDF(fileinfo.Filename,product_name,[0 1],[-9 1]);

% 基于Matlab自带函数读取
% lat_s = hdfread(fileinfo.Filename, '/Latitude');
% lat_s = hdfread(fileinfo.Filename, '/Latitude', 'Index', {[1  2],[1  1],[3744 1]});

% 在指定纬度范围内确定廓线的索引
latRangeWanted = lat > latMin & lat < latMax;
firstProfile = find(latRangeWanted,1,'first');
lastProfile = find(latRangeWanted,1,'last');
numberProfiles = lastProfile - firstProfile + 1;

%% 读取消光系数及其不确定度
product_name = 'Extinction_Coefficient_532';
[info, sigma] = readHDF(fileinfo.Filename,product_name,[firstProfile-1 0],[numberProfiles -9]);

% 基于Matlab自带函数读取
% sigma_s = hdfread(fileinfo.Filename, '/Extinction_Coefficient_532');
% sigma_s = hdfread('CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf', '/Extinction_Coefficient_532', 'Index', {[firstProfile 1],[1  1],[numberProfiles 399]});

product_name = 'Extinction_Coefficient_Uncertainty_532';
[info, sigma_uncer] = readHDF(fileinfo.Filename,product_name,[firstProfile-1 0],[numberProfiles -9]);

% 基于Matlab自带函数读取
% sigma_uncer_s = hdfread(fileinfo.Filename, '/Extinction_Coefficient_Uncertainty_532');
% sigma_uncer_s = hdfread('CAL_LID_L2_05kmAPro-Prov-V3-01.2010-04-17T20-09-13ZN.hdf', '/Extinction_Coefficient_Uncertainty_532', 'Index', {[firstProfile 1],[1  1],[numberProfiles 399]});

% 矩阵转置，将sigma与sigma_uncer统一成 numberAlts x numberProfiles
sigma = sigma';
sigma_uncer = sigma_uncer';

%% 读取质量评价数据
% QC, CAD_Score, AVD 是三维矩阵,如：2 x numberAlts x numberProfiles
product_name = 'Extinction_QC_Flag_532';
[info, extQC] = readHDF(fileinfo.Filename,product_name,[firstProfile-1 0 0],[numberProfiles -9 -9]);  

% extQC_s = hdfread(fileinfo.Filename, 'Extinction_QC_Flag_532');
% extQC_s = hdfread(fileinfo.Filename, 'Extinction_QC_Flag_532', 'Index', {[firstProfile 1 1],[1 1 1],[numberProfiles 399 2]});
% extQC_s = permute(extQC_s,[3 2 1]);

product_name = 'CAD_Score';
[info, CAD] = readHDF(fileinfo.Filename,product_name,[firstProfile-1 0 0],[numberProfiles -9 -9]);     

% CAD_s = hdfread(fileinfo.Filename,'/CAD_Score');
% CAD_s = hdfread(fileinfo.Filename,'/CAD_Score', 'Index', {[firstProfile 1 1],[1 1 1],[numberProfiles 399 2]});
% CAD_s = permute(CAD_s,[3 2 1]);

product_name = 'Atmospheric_Volume_Description';
[info, AVD] = readHDF(fileinfo.Filename,product_name,[firstProfile-1 0 0],[numberProfiles -9 -9]); 

% AVD_s = hdfread(fileinfo.Filename,'/Atmospheric_Volume_Description');
% AVD_s = hdfread(fileinfo.Filename,'/Atmospheric_Volume_Description', 'Index', {[firstProfile 1 1],[1 1 1],[numberProfiles 399 2]});
% AVD_s = permute(AVD_s,[3 2 1]);

%% 筛选指定数据
% 使用逻辑数组，筛选出指定数据。保留均值为“真”的元素，丢弃“错误”的元素，将索引结果存储在矩阵useSamples中。
useSamples = true(numberAlts,numberProfiles);

% 筛选掉填充值
useSamples(sigma == -9999) = false;

%% 由于星载激光雷达的信噪比低，需要对星载激光雷达观测的气溶胶进行筛选
% 基于气溶胶特征，筛选非气溶胶数据
[ftype, subtype] = get_feature_type_and_subtype(AVD);

% AVD的大小为2 x numberAlts x numberProfiles，而sigma的大小为numberAlts x numberProfiles。
% 这是因为CALIPSO在平流层记录的消光系数垂直分辨率是60米，而在对流层记录的是30米分辨率，所以采用了三维矩阵。
% 参阅CALIPSO Data Qualtity Summaries中的Atmospheric Volume Description条目。
% 检查30米分辨率和60米分辨率包含的气溶胶特征，如果包含则逻辑值为真，否则为假。对于气溶胶，特征类型 = 3。

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

%% 基于CAD Score，筛选非气溶胶数据
% 对于气溶胶类型的分类，CAD反映了反演的气溶胶可信度，只有在CAD>-20的情况下，产品分类的可信度才高。
% 这里将筛选出CAD score 大于-20的气溶胶层。
% For this example, we will screen out aerosol layers whose CAD score is
% greater than -20. For a conservative approach, we will screen out the
% entire 60 m sigma bin if either 30 m element in that bin has CAD > -20.
% Furthermore, since cloud CAD scores are also reported in the profile
% product CAD score arrays, we will first make all cloud scores equal to
% fill values so we don't screen out cases where a 30 m cloud and a 30 m
% aerosol layer (with CAD < -20) exist in the same 60 m element.

CAD_fill = -127;
CAD(CAD > 0) = CAD_fill;
goodCAD = squeeze( (CAD(1,:,:) < -20 | CAD(1,:,:) == CAD_fill) &...
                   (CAD(2,:,:) < -20 | CAD(2,:,:) == CAD_fill) );
useSamples(~goodCAD) = false;

%% 基于QC flags，筛选非气溶胶数据
% 只选用反演过程中初始雷达比在反演迭代过程未变（QC = 0）的情况，
% 以及激光雷达探测到的气溶胶上方和下方大气分子散射强度能得到的情况（QC = 1）。

% The CALIPSO extinction QC flags summarize the initial and final state of
% the extinction retrieval. During the iterative process to find a solution
% the lidar ratio may be adjusted. Often, we have most confidence in 
% extinction solutions when the lidar ratio is unchanged during the 
% retrieval (extinction QC = 0) or if the retrieval is constrainted
% (extinction QC = 1). See the final lidar ratio and extinction QC entries
% in the CALIPSO Layer Products Data Quality Summaries for more details.

% The extinction QC flags reported are reported from both clouds and
% aerosols in the array. Because the extinction is retrieved from both 30 m
% bins in the 60 m bin, regardless of feature type, we will check the
% extinction QC flags for all features (including clouds in the array). If
% either 30 m element is not 0 or 1, we will throw out the whole 60 m bin.
extQCfill = 32768;
extQC(ftype ~= 3) = extQCfill;
goodExtQC = squeeze(...
    (extQC(1,:,:) == 0 | extQC(1,:,:) == 1 | extQC(1,:,:) == extQCfill) & ...
    (extQC(2,:,:) == 0 | extQC(2,:,:) == 1 | extQC(2,:,:) == extQCfill)); 
useSamples(~goodExtQC) = false;

%% 基于extinction uncertainty，筛选非气溶胶数据
% CALIPSO消光不确定性是在迭代过程中求解，有时激光雷达方程无法收敛。
% 当检索到的消光开始无限增长时，消光不确定性被赋值为99.99。
% 当消光不确定度为99.99时，消光值不应被信任。
useSamples(sigma_uncer > 99.9) = false;

%% 绘制消光剖面

meanSigmaUnscreened = zeros(numberAlts,1);
meanSigmaUnscreened_uncer = zeros(numberAlts,1);
numSamplesUnscreened = zeros(numberAlts,1);
meanSigmaScreened = zeros(numberAlts,1);
meanSigmaScreened_uncer = zeros(numberAlts,1);
numSamplesScreened = zeros(numberAlts,1);

for k = 1:numberAlts

    t = sigma(k,:) ~= -9999;
    meanSigmaUnscreened(k) = mean(sigma(k,t));
    meanSigmaScreened(k) = mean(sigma(k,useSamples(k,:)));

    % 通过传播消光不确定性来计算不确定性，即不确定性的平方是误差的平方和除以样本数量的总和。
    numSamplesUnscreened(k) = sum(t);
    if numSamplesUnscreened(k) > 0
        meanSigmaUnscreened_uncer(k) = sqrt(sum(sigma_uncer(k,t).^2)./numSamplesUnscreened(k));
    end
    numSamplesScreened(k) = sum(useSamples(k,:));
    if numSamplesScreened(k) > 0
        meanSigmaScreened_uncer(k) = sqrt(sum(sigma_uncer(k,useSamples(k,:)).^2)./numSamplesScreened(k));
    end
end

% -- Plot of mean extinction vs. altitude, unscreened and screened

figure('position',[1023 300 464 681],'paperpositionmode','auto')
plot(meanSigmaUnscreened,alts,'linewidth',2,'linestyle','--')
hold on
plot(meanSigmaScreened,alts,'linewidth',2,'color','red')
set(gca,'ylim',[-0.5 8])
ylabel('Altitude (km)')
xlabel('Mean sigma (km^{-1})')
title('CALIPSO Mean Extinction 532 nm')
legend('Unscreened','Screened')
grid on

% -- Plot of relative uncertainty vs. altitude with uncertainty, unscreened and screened

figure('position',[1023 300 464 681],'paperpositionmode','auto')
relUncertUnscreened = 100*(meanSigmaUnscreened_uncer./meanSigmaUnscreened);
relUncertScreened = 100*(meanSigmaScreened_uncer./meanSigmaScreened);
plot(relUncertUnscreened,alts,'linewidth',2,'linestyle','--')
hold on
plot(relUncertScreened,alts,'linewidth',2,'color','red')
set(gca,'ylim',[-0.5 8])
set(gca,'xlim',[50 400])
ylabel('Altitude (km)')
xlabel('Relative extinction uncertainty (%)')
title('CALIPSO Extinction 532 nm Relative Uncertainty')
legend('Unscreened','Screened')
grid on

% -- Plot of number of samples vs. altitude, unscreened and screened

figure('position',[1023 300 464 681],'paperpositionmode','auto')
plot(numSamplesUnscreened,alts,'linewidth',2,'linestyle','--')
hold on
plot(numSamplesScreened,alts,'linewidth',2,'color','red')
set(gca,'ylim',[-0.5 8])
ylabel('Altitude (km)')
xlabel('Number')
title('CALIPSO Number of Samples')
legend('Unscreened','Screened')
grid on