ENVI IDL：如何生成FY4A快照

news2025/1/11 22:53:22

01 数据说明

在这里插入图片描述

FY4A全圆盘（DISK，全球）多光谱影像，panoply软件打开数据层次结构如下：

在这里插入图片描述

我们生成快照主要使用到其中的NOMChannel01、NOMChannel02、NOMChannel03进行快照显示，注意我并没有进行辐射定标。

02 生成JPG快照

代码如下：

; @Author	: ChaoQiezi
; @Time		: 2023年11月13日-上午8:40:16
; @Email	: chaoqiezi.one@qq.com

; 该程序用于 生成风云4A快照(.png)

function read_band, h5_path, ds_name, range_name
    ds = read_h5(h5_path, ds_name=ds_name, /double)
    ds_range = read_h5(h5_path, ds_name=ds_name, attr_name=range_name, /double)
    invalid_pos = where((ds lt ds_range[0]) or (ds gt ds_range[1]), /null)
    ds[invalid_pos] = !values.F_NAN
    
    return, ds
end

pro fy4a_quicklook
    ; 准备
    in_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\Week8\FY4A\'
    out_dir = 'D:\Objects\JuniorFallTerm\IDLProgram\Experiments\ExperimentalData\Week8\out_fy4a_quicklook\'
    if ~file_test(out_dir, /directory) then file_mkdir, out_dir
    band1_name = 'NOMChannel01'
    band2_name = 'NOMChannel02'
    band3_name = 'NOMChannel03'
    range_name = 'valid_range'
    
    ; 检索和循环
    h5_paths = file_search(in_dir, '*DISK*.hdf')
    foreach h5_path, h5_paths do begin
        ; 数据集和属性
        band1 = read_band(h5_path, band1_name, range_name)
        band2 = read_band(h5_path, band2_name, range_name)
        band3 = read_band(h5_path, band3_name, range_name)
        bands = make_array([3, size(band1, /dimensions)], /uint)
        
        band1 = bytscl(band1, /nan)
        band2 = bytscl(band2, /nan)
        band3 = bytscl(band3, /nan)
        
        bands[0, *, *] = band3
        bands[1, *, *] = band2
        bands[2, *, *] = band1
        
        ; 输出png
        jpg_path = out_dir + file_basename(h5_path, '.hdf') + '.jpg'
        png_path = out_dir + file_basename(h5_path, '.hdf') + '.png'
        write_jpeg, jpg_path, bands, true=1, /order
    endforeach
end

里面涉及的read_h5函数我之前一直在封装，所以移植不太方便了，但是还是贴出来，大家自己写一些比较简单的会更好：

;+
;   函数用途:
;       (私有函数)用于读取HDF5文件指定ID(HDF文件ID或数据集ID)下的属性
;   函数参数:
;       h5_id: HDF5文件ID或数据集ID
;       attr_name: 属性名称
;-
function _read_h5_attribute, h5_id, attr_name
    ; 获取属性ID
    attr_id = h5a_open_name(h5_id, attr_name)
    
    ; 获取属性
    attr = h5a_read(attr_id)
    h5a_close, attr_id
    
    return, attr
end

;+
;   函数用途:
;       (私有函数)用于读取HDF5文件下的数据集
;   函数参数:
;       h5_id: HDF5文件ID
;       ds_name: 数据集名称
;-
function _read_h5_dataset, h5_id, ds_name
    ; 获取数据集
    ds_id = h5d_open(h5_id, ds_name)
    data = h5d_read(ds_id)
    h5d_close, ds_id
    
    return, data
end

;+
;   函数用途:
;       该函数用于读取HDF5文件的数据集及其属性
;   函数参数:
;       h5_path: HDF5文件的输入路径
;       ds_name(关键字传参): 读取数据集的路径(注意: 路径为文件内路径)
;       attr_name(关键字传参): 读取属性的名称
;       double(关键字传参): 输出数据的类型是否为双精度浮点
;-
function read_h5, h5_path, ds_name=ds_name, attr_name=attr_name, double=double
    ; 获取HDF5文件数据集ID
    h5_id = h5f_open(h5_path)
    
    if keyword_set(ds_name) && keyword_set(attr_name) then begin
        ds_id = h5d_open(h5_id, ds_name)
        data = _read_h5_attribute(ds_id, attr_name)
        h5d_close, ds_id
    endif else if keyword_set(ds_name) then begin
        data = _read_h5_dataset(h5_id, ds_name)
    endif else if keyword_set(attr_name) then begin
        data = _read_h5_attribute(h5_id, attr_name)
    endif else begin
        return, -1
    endelse
    
    if keyword_set(double) then data = double(data)
    
    ; 关闭资源
    h5f_close, h5_id
    
    return, data
end

封装的比较复杂了，需要自己多多体会了。

里面涉及的关键只有两个，一个是write_jpeg函数，一个是bytscl函数。

2.1 `write_jpeg`函数说明：

在这里插入图片描述

/order参数表示绘制的顺序，默认就是从上到下绘制。如果设置了就是从下往上绘制，我们的栅格矩阵一般最上面就是北边，正常的，不需要设置该参数。
true表示输入Image数组(三维)中波段维度的位置。如果Image数组的形式为(波段数，列数，行数)，那么true为1；如果Image数组的形式为(列数，波段数，行数)，那么true为2；如果Image数组的形式为(列数，行数，波段数)，那么true为3；
quality表示压缩的质量，默认是75的压缩。