1. 实验内容
这里没什么新的知识点,就是各个地方的数组定义要注意一下啊,有的是无符号数组,有的是字节型数组。最变态的是由于风云四号卫星数据比较大,一个数组占用内存甚至达到300M+,所以不用的内存需要即使处理,否则如果是循环那么很容易出现内存无法分配到数组的错误出现。还有就是一些函数你可能了解,但是有些参数陌生所以需要多看帮助了解了解。
总体思路:
1. 打开目录获取相关数据。
2. 对获取的数据进行提取获取三个波段的数据集。
3. 然后使用三维数组存储上面的三个波段数据并用write_tiff输出即可。
上面是波段合成,下面是线性拉伸。
1. 将范围限制在2%~98%之间。
2. 将一些细节进行处理,如南北极转换、数组类型创建正确等等。
3. 借write_jpeg函数输出即可。
2.编程
2.1 代码部分
function get_hdf5_ds, file_path, ds_name ; 该函数是获取hdf4文件的数据集 ; 获取文件id file_id = h5f_open(file_path) ; 获取数据集id ds_id = h5d_open(file_id, ds_name) ; 获取数据集数据 ds_data = h5d_read(ds_id) ; 关闭文件 h5f_close, file_id ; 返回数据集数据 return, ds_data end pro week_seven_test ; 本程序用于波段合成和线性拉伸(以风云4号数据为例) ; 这里需要说一嘴,风云4号卫星数据我现在拿的是hdf5数据,但是它的后缀是.hdf而不是.he5,搞得我前面的函数写错了(现在订正过来了) ; 路径 in_dir = 'D:/IDL_program/experiment_data/chapter_4/' out_dir = 'D:/IDL_program/experiment_data/chapter_4/tiff_520/' if file_test(out_dir, /directory) eq 0.0 then file_mkdir, out_dir ; 不存在该目录就创建 ; 获取in_dir路径下的hdf5文件的所有路径 hdf_path_list = file_search(in_dir, '*DISK*.hdf', count=hdf_count) ; 注意这个hdf5文件不是he5作为后缀名 ; 获取hdf5文件的所有名称 hdf_name_list = file_basename(hdf_path_list, '.hdf') ; 去掉了后缀名 ; 循环处理hdf文件 for hdf_i = 0, hdf_count - 1 do begin ; 获取该循环下的hdf文件的路径 hdf_path = hdf_path_list[hdf_i] ; 获取该循环下的hdf文件的名称 hdf_name = hdf_name_list[hdf_i] ; 输出路径 out_path = out_dir + hdf_name + '_mix.tiff' ; 获取hdf文件的三个数据集的数据(由hdf explorer知道该数据是16-bit-unsigned-integer) band_red = get_hdf5_ds(hdf_path, 'NOMChannel03') band_green = get_hdf5_ds(hdf_path, 'NOMChannel02') band_blue = get_hdf5_ds(hdf_path, 'NOMChannel01') ; 获取该数据的size band_size = size(band_red) ; col band_col = band_size[1] ; row band_row = band_size[2] ; 创建合成波段的数组 band = uintarr(3, band_col, band_row) ; 你可能会uintarr(band_col, band_row, 3),但是后面的函数要求是要波段数放在前面 ; 将三个波段放在上面的波段合成数组中 band[0, *, *] = band_red band[1, *, *] = band_green band[2, *, *] = band_blue ; 后面都不会使用这几个变量了,但是它们占据的内存空间又特别大,所以将它们赋空值释放空间 band_red = !null band_green = !null band_blue = !null ; 将无效值(有效值是0-4095)变成0 band = ((band lt 4095) and (band gt 0)) * band ; 波段合成并输出 write_tiff, out_path, band, /short ; 你将这个band传入后write_tiff函数会自己进行波段合成 ; 另外,由于band数组是uintarr,所以指定数据类型是/short ; 线性拉伸 ; 获取一个波段的数组元素个数 band_element = n_elements(band[0, *, *]) ; 知道一个波段的数组2% 98%的元素个数 min_pos = 0.02 * band_element max_pos = 0.98 * band_element ; 将每一个波段的数据进行排序 ; 创建存储拉伸后的数据的数组 jpeg_data = bytarr(3, band_col, band_row) ; 需要字节型数组 for band_i = 0, 2 do begin ; 获取该循环下的波段数据 temp_band = band[band_i, *, *] ; sort()将排序(从小到大)之后的索引按顺序返回 band_sort = sort(temp_band) ; 获取最大最小值的索引 min_index = band_sort[min_pos] max_index = band_sort[max_pos] ; 获取最大最小值 min_data = temp_band[min_index] max_data = temp_band[max_index] ; 对其进行拉伸 jpeg_data[band_i, *, *] = bytscl(temp_band, min=min_data, max=max_data) ; 解除内存占用 temp_band = !null endfor ; 输出为jpeg图像文件 write_jpeg, out_dir + hdf_name_list[hdf_i] + '_way.jpeg', jpeg_data,true=1, order=1 ; true=1表示传入数据中波段数在第一维,2表示在第二维,3表示在第三维。order表示翻转图片(由于做完之后会南北极颠倒所以需要翻转过来) ; 关闭内存占用 jpeg_data = !null band = !null endfor print, '>>> 完毕!' end
2.2 运行结果
黑图是没有经过线性拉伸的,所以普通照片查看器只能隐隐看得到一些(但是使用ENVI、Arcgis可以清晰查看)。清晰的地球照片是经过线性拉伸的。可以直接用照片查看器打开