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的压缩。
2.2 bytscl
函数说明:
这个函数非常关键,由于我们自己的数组可能是浮点型或者双精度等,在生成JPEG
时实际上会导致输出的图片特别失真或者直接全黑或者全白等,我们需要将数组转换为字节型数组,即范围是0~255的整型范围。
bytscl
函数通过下方公式进行类似归一化的缩放到0~255:
默认Top
是255。
/nan
在数组中存在NAN无效值时使用。
最终处理结果如下:
仅仅用来查看一下FY4A,不可用于数据分析,因为像元值完全发生变化,经过映射到0~255等操作。