前言
这是有关充分利用地理空间应用程序的技巧系列的一部分。使用IDL 应用程序打开以 HDF5 格式存储的 HICO 数据集。该文章同样描述了如何在 ENVI 中使用 H5_Browser 或新的 HDF5 Reader 打开 HDF5 格式的 HICO 文件。
本博客演示了如何实现 IDL 代码以将澳大利亚蒙哥马利礁的 HDF5 HICO 场景打开为 ENVI 格式。 随后的步骤包括为进一步分析准备结果数据。
1、代码准备
本例中使用的 HICO 数据集 (H2012095004112.L1B_ISS) 是从 NASA GSFC 档案下载的,可以通过俄勒冈州立大学的 HICO 网站或 NASA GSFC Ocean Color 网站访问。 请注意,您还可以通过 OSU 网站申请成为注册的 HICO 数据用户,从而获得对 ENVI 格式数据集的访问权限。
本示例中使用的 IDL 代码可从 NASA GSFC Ocean Color 网站的 Documents > Software/Tools > IDL Library > hico 获得。 您需要的三个 IDL 文件是:byte_ordering.pro、nrl_hico_h5_to_flat.pro 和 write_nrl_header.pro。但上述代码的下载链接已经失效,我把它的原版找到,放到下面
nrl_hico_h5_to_flat
pro nrl_hico_h5_to_flat,infile,outdir
compile_opt idl2
program_name = 'nrl_hico_h5_to_flat'
program_version = '1.1'
program_date = '2013/11/05'
program_author = 'Marcos Montes, marcos.montes@nrl.navy.mil'
if (n_params() ne 2) then begin
if (n_elements(infile) ne 1) then begin
infile=!null
infile=dialog_pickfile(/must_exist,filter='*.h5',title='Select a HICO HDF5 file')
if (infile eq '') then begin
print,'==============================================='
print,'nrl_hico_h5_to_flat: input HDF5 file not selected; returning to caller'
print,'==============================================='
return
endif
endif
if (n_elements(outdir) ne 1) then begin
outdir=!null
path=file_dirname(infile,/mark_directory)
outdir=dialog_pickfile(path=path,/directory,/must_exist,title='Select an output directory')
if (outdir eq '') then begin
print,'==============================================='
print,'nrl_hico_h5_to_flat: output directory not selected; returning to caller'
print,'==============================================='
return
endif
endif
; print,'==============================================='
; print,'In : ',program_name
; print,'Returning because number of provided arguments NE 2'
; print,''
; print,'Usage:'
; print,'IDL> nrl_hico_h5_to_flat,infile,outdir'
; print,'2 parameters required. First is input h5 file'
; print,'second is an output directory'
; print,'==============================================='
; return
endif
if (~h5f_is_hdf5(infile)) then begin
print,'==============================================='
print,'In : ',program_name
print,'Returning; infile is NOT an HDF5 file, infile = ',infile,$
format='(2a)'
print,'Please ensure that the directory path and name are correct'
print,'==============================================='
return
endif
if (~file_test(outdir,/directory,/write)) then begin
print,'==============================================='
print,'Returning; outfir is not a writeable directory, outdir = ',outdir,$
format='(2a)'
print,'==============================================='
return
endif
; This IS overkill. Basically, I'm going to read the entire file at
; once, and import all the data into structures. Hey, it works on my
; MacBook Pro with my version of ENVI, so it will work for everyone, right?
; The benefit is that everything is done at once. The (possible)
; disadvantage is that everything is done at once. I will endeavor to test
; this on a few different OS/IDL combinations to see if it is robust
; enough. Thankfully, HICO files are a standard size, which puts a limit on
; the amount of data we expect.
hparse=h5_parse(infile,/read_data)
base=file_basename(hparse._name,'.h5')
this_bo=byte_ordering()
; BEGIN RADIANCE FILE
; First, the calibrated HICO data and associated header file
dims=hparse.products._lt._Dimensions
oname=string(outdir,path_sep(),base,'.bil',format='(4a)')
openw,lun,oname,/get_lun
writeu,lun,transpose(hparse.products._lt._data,[1,0,2]) ; re-ordered to BIL from BIP
close,lun & free_lun,lun
isf=1.d/reform(hparse.products._lt.slope._data)
isf_array=make_array(dims[0],type=4,value=floor(isf)*1.0) ; scale factors
runits=reform(hparse.products._lt.units._data) ;units
zplot_y=string(isf,'* radiance (',runits,')',format='(f5.2,3a)') ; zplot
scene_id=reform(hparse.metadata.hico.target.id._data)
target_name=reform(hparse.metadata.hico.target.name._data)
target_type=reform(hparse.metadata.hico.target.type._data)
target_aos=reform(hparse.metadata.hico.target.aos._data)
target_angle=reform(hparse.metadata.hico.target.angle._data)
description=['Calibrated HICO data',$
string(' target_name: ',target_name),$
string(' target_type: ',target_type),$
string(' target_angle: ',target_angle),$
string(' target_aos: ',target_aos)]
; construction of history
history=['']
tnames=tag_names(hparse.metadata.hico.calibration)
lens=strlen(tnames)
for i=6,68 do begin ; normal radiance history part
key=strlowcase(strmid(reform(tnames[i]),0,lens[i]-1))
value=strtrim(reform(hparse.metadata.hico.calibration.(i).(3)),2)
history=[history,string(key,' = ',value,format='(3a)')]
endfor
history=[history,'']
tnames=tag_names(hparse.metadata.hico.target)
lens=strlen(tnames)
for i=6,10 do begin ; target_* part
key=string('target_',strlowcase(strmid(reform(tnames[i]),0,lens[i])),format='(2a)')
value=strtrim(reform(hparse.metadata.hico.target.(i).(3)),2)
history=[history,string(key,' = ',value,format='(3a)')]
endfor
history=[history,'']
tnames=tag_names(hparse.metadata.fgdc.identification_information.time_period_of_content)
lens=strlen(tnames)
for i=6,9 do begin ; "date/time" part
key=string(strlowcase(strmid(reform(tnames[i]),0,lens[i])),format='(a)')
value=strtrim(reform(hparse.metadata.fgdc.identification_information.time_period_of_content.(i).(3)),2)
history=[history,string(key,' = ',value,format='(3a)')]
endfor
history=[history,'']
tnames=tag_names(hparse.metadata.fgdc.identification_information.spatial_domain)
lens=strlen(tnames)
for i=6,9 do begin ; "sort of like bounding box part"
key=string(strlowcase(strmid(reform(tnames[i]),0,lens[i])),format='(a)')
value=strtrim(reform(hparse.metadata.fgdc.identification_information.spatial_domain.(i).(3)),2)
history=[history,string(key,' = ',value,format='(3a)')]
endfor
history=[history,'']
history=[history,string('hdf_extraction_program_name_and_version = ',program_name,' ',program_version,$
format='(4a)')]
history=[history,string('hdf_extraction_program_author = ',program_author,format='(2a)')]
history=[history,string('input_hdf_file = ',infile,format='(2a)')]
history=[history,'']
oname=string(outdir,path_sep(),base,'.hdr',format='(4a)')
openw,lun,oname,/get_lun
write_nrl_header,lun,ns=dims[1],nl=dims[2],nb=dims[0],offset=0,inter='bil',$
data_type=2,file_type='ENVI Standard',sensor_type='HICO',$
wl=0.001d*hparse.products._lt.wavelengths._data,$
wavelength_unit='micrometers',$
fwhm=0.001d*hparse.products._lt.fwhm._data,$
units=runits,xstart=1,ystart=1,$
scale=isf_array,byte_order=this_bo,$
z_titles=['wavelength (micrometer)',zplot_y],$
scene_id=scene_id,descrip=description,$
history=history
close,lun & free_lun,lun
; Note: currently the header file is missing some of the values. It seems
; to me that the L0 Hico Header data isn't present; but it is more
; likely that I haven't yet figured out exactly where the data is. I
; see the keywords; I don't see the values.
; END RADIANCE FILE
; BEGIN GEOMETRY FILE
; Now a geometry file; convert from 6 BSQ planes to BIL
oname=strtrim(reform(hparse.metadata.hico.inputs.geometry_file.(3)),2)
oname=string(outdir,path_sep(),oname,format='(3a)')
nlen=strlen(oname)
odat=transpose([[[hparse.navigation.(7).(5)]],[[hparse.navigation.(6).(5)]],$
[[hparse.navigation.(9).(5)]],[[hparse.navigation.(8).(5)]],$
[[hparse.navigation.(11).(5)]],[[hparse.navigation.(10).(5)]]],[0,2,1])
openw,lun,oname,/get_lun
writeu,lun,odat
close,lun & free_lun,lun
odat=!null ; free the space
; order is such to get the bands listed below:
band_names=['longitude (degrees)', 'latitude(degrees)', $
'view zenith (degrees)','view azimuth (degrees clockwise from N)', $
'solar zenith (degrees)', 'solar azimuth (degrees clockwise from N)']
oname=string(strmid(oname,0,nlen-3),'hdr',format='(2a)')
openw,lun,oname,/get_lun
write_nrl_header,lun,ns=dims[1],nl=dims[2],nb=6,offset=0,inter='bil',$
data_type=4,file_type='ENVI Standard',xstart=1,ystart=1,$
byte_order=this_bo,bnames=band_names,sensor_type='HICO',$
history=history
close,lun & free_lun,lun
; END GEOMETRY FILE
end
byte_ordering
function byte_ordering
; Want to determine byte ordering on this machine.
; Not sure how to this in IDL, so I will use the same code as Fortran
; Returns ENVI's sense of byte ordering
; No inputs
; usage:
; IDL> a=byte_ordering()
compile_opt idl2
openw,lun,'junk0912875634.bin',/get_lun,/delete
one=1s ; must be 2 byte integer
hold=bytarr(2) & hold[*]=byte(0)
writeu,lun,one ; write a 2-byte integer
point_lun,lun,0 ; rewind
readu,lun,hold ; read as 2@1byte values
close,lun & free_lun,lun ; cleanup
return,fix(hold[1]) ; return function's value
end
write_nrl_header
Pro write_nrl_header, out_unit, ns=ns, nl=nl, nb=nb, offset=offset, $
inter=inter, data_type=data_type, $
byte_order=byte_order, bbl=bad_bands, $
file_type=file_type, sensor_type=sensor_type, $
bnames=band_names, descrip=descrip, wl=wave_len, $
fwhm=wave_fwhm, $
wavelength_unit=wavelength_unit,$
xstart=xstart, ystart=ystart, map_info=map_info, $
def_bands=def_bands, $
z_range=z_range, z_average=z_average, $
z_titles=z_titles, $
def_stretch=def_stretch, $
pixel_size=pixel_size, $
scale=image_scale_factor, $
units=image_unscaled_units, $
date=image_center_date, $
time=image_center_time, $
longitude=image_center_long, $
long_hem=image_center_long_hem, $
latitude=image_center_lat, $
lat_hem=image_center_lat_hem, $
zenith_ang=image_center_zenith_ang, $
azimuth_ang=image_center_azimuth_ang, $
altitude=sensor_altitude, $
tafkaa_ground_elevation=ground_elevation,$
image_data_quantity=image_data_quantity,$
history=history,noenvi=noenvi,$
version=version,scene_id=scene_id,$
missing_lines=missing_lines,$
num_missing_lines=num_missing_lines,$
spectra_names=spectra_names,$
reflectance_scale_factor=reflectance_scale_factor,$
major_frame_offsets=major_frame_offsets,$
minor_frame_offsets=minor_frame_offsets
if (keyword_set(noenvi)) then begin
endif else printf, out_unit, format='("ENVI")'
if n_elements(descrip) ne 0 then begin
printf,out_unit,'description = { ',descrip,format='(a,(a," "),$)'
printf,out_unit,'}'
endif
if n_elements(ns) ne 0 then printf, out_unit, format='("samples = ", i0)', ns
if n_elements(nl) ne 0 then printf, out_unit, format='("lines = ", i0)', nl
if n_elements(nb) ne 0 then printf, out_unit, format='("bands = ", i0)', nb
if n_elements(offset) ne 0 then printf, out_unit, format='("header offset = ", i0)', offset
if n_elements(inter) ne 0 then printf, out_unit, format='("interleave = ", a)', inter
if n_elements(data_type) ne 0 then printf, out_unit, format='("data type = ", i0)', $
data_type
if n_elements(byte_order) ne 0 then printf, out_unit, format='("byte order = ", i0)', $
byte_order
if n_elements(file_type) ne 0 then printf, out_unit, format='("file type = ", a)', $
file_type
if n_elements(sensor_type) ne 0 then printf, out_unit, format='("sensor type = ", a)', $
sensor_type
if n_elements(version) ne 0 then printf, out_unit, format='("version = ", i0)', $
version
if n_elements(scene_id) ne 0 then printf, out_unit, format='("scene_id = ", i0)', $
scene_id
if n_elements(xstart) ne 0 then printf, out_unit, format='("x start = ", i0)', xstart
if n_elements(ystart) ne 0 then printf, out_unit, format='("y start = ", i0)', ystart
if (n_elements(reflectance_scale_factor) ne 0) then begin ; MJM 2010/09/22
printf,out_unit,format='("reflectance scale factor = ", f0)',$
reflectance_scale_factor
endif
if (n_elements(minor_frame_offsets) ne 0) then begin ; MJM 2010/09/22
printf,out_unit,format='(a,i0,",",i0,a)','minor frame offsets = { ',$
minor_frame_offsets,'}'
endif
if (n_elements(major_frame_offsets) ne 0) then begin ; MJM 2010/09/22
printf,out_unit,format='(a,i0,",",i0,a)','major frame offsets = { ',$
major_frame_offsets,'}'
endif
if (n_elements(num_missing_lines) ne 0) then begin ; MJM 2010/09/22
printf,out_unit,format='("num_missing_lines = ",i0)',num_missing_lines
endif
if (n_elements(missing_lines) ne 0) then begin ; MJM 2010/Sep/22, modeled on bbl
printf, out_unit, format='(a)','missing_lines = {'
printf, out_unit, format='(10(" ",i0,:,","),$)', missing_lines
printf, out_unit, format='("}")'
endif
if n_elements(map_info) ne 0 then begin
num = n_elements(map_info)
text = '{'
for i = 0, num-2 do begin
text = text + map_info[i] + ', '
endfor
text = text + map_info[num-1] + '}'
printf, out_unit, format='("map info = ",a)', text
endif
if n_elements(def_bands) ne 0 then begin
printf, out_unit, format='("default bands = {",3(i0,:,", "),$)', def_bands
printf, out_unit, format='("}")'
endif
if n_elements(z_range) ne 0 then $
printf, out_unit, format='("z plot range = {",f0.2,", ",f0.2,"}")', z_range
if n_elements(z_average) ne 0 then $
printf, out_unit, format='("z plot average = {",i0,", ",i0,"}")', z_average
if n_elements(z_titles) ne 0 then $
printf, out_unit, format='("z plot titles = {",a,", ",a,"}")', z_titles
if n_elements(def_stretch) ne 0 then $
printf, out_unit, format='("default stretch = ", a)', $
def_stretch
if n_elements(pixel_size) ne 0 then begin
num = n_elements(pixel_size)
text = '{'
for i = 0, num-2 do begin
text = text + pixel_size[i] + ', '
endfor
text = text + pixel_size[num-1] + '}'
printf, out_unit, format='("pixel size = ",a)', text
endif
if n_elements(image_unscaled_units) ne 0 then begin
num = n_elements(image_unscaled_units)
text = '{'
for i = 0, num-2 do begin
text = text + image_unscaled_units[i] + ', '
endfor
text = text + image_unscaled_units[num-1] + '}'
printf, out_unit, format='("image_unscaled_units = ",a)', text
endif
if (n_elements(wavelength_unit) ne 0) then $
printf,out_unit, format='("wavelength units = ",a)',wavelength_unit
if (n_elements(image_data_quantity) ne 0) then $
printf, out_unit, 'image_data_quantity = ',image_data_quantity,format='(2a)'
if n_elements(image_center_date) ne 0 then $
printf, out_unit, format='("image_center_date = {",i4.4,", ",i2.2,", ",i2.2,"}")', $
image_center_date
if n_elements(image_center_time) ne 0 then $
printf, out_unit, format='("image_center_time = {",i2.2,", ",i2.2,", ",f8.3,"}")', $
image_center_time
if n_elements(image_center_long) ne 0 then $
printf, out_unit, format='("image_center_long = {",f0,", ",f0,", ",f8.3,"}")', $
image_center_long
if n_elements(image_center_long_hem) ne 0 then $
printf, out_unit, format='("image_center_long_hem = ", a)', image_center_long_hem
if n_elements(image_center_lat) ne 0 then $
printf, out_unit, format='("image_center_lat = {",f0,", ",f0,", ",f8.3,"}")', $
image_center_lat
if n_elements(image_center_lat_hem) ne 0 then $
printf, out_unit, format='("image_center_lat_hem = ", a)', image_center_lat_hem
if n_elements(image_center_zenith_ang) ne 0 then begin
printf, out_unit, format='("image_center_zenith_ang = {",f0,", ",f0,", ",f8.3,"}")', $
image_center_zenith_ang
endif
if n_elements(image_center_azimuth_ang) ne 0 then begin
printf, out_unit, format='("image_center_azimuth_ang = {",f0,", ",f0,", ",f8.3,"}")', $
image_center_azimuth_ang
endif
if n_elements(ground_elevation) ne 0 then begin
printf,out_unit,format='("tafkaa_ground_elevation = ",f0)',ground_elevation
endif
if n_elements(sensor_altitude) ne 0 then begin
printf, out_unit, format='("sensor_altitude = ",f9.2)', sensor_altitude
endif
if n_elements(image_scale_factor) ne 0 then begin
printf, out_unit, format='("image_scale_factor = {")'
printf, out_unit, format='(7(" ",f0,:,","),$)', image_scale_factor
printf, out_unit, format='("}")'
endif
if n_elements(band_names) ne 0 then begin
printf, out_unit, format='("band names = {")'
printf, out_unit, format='(7(" ",a,:,","),$)', band_names
printf, out_unit, format='("}")'
endif
if (n_elements(spectra_names) ne 0) then begin ; MJM 2010/Sep/22, modeled on band names, but 3 per line
printf, out_unit, format='("spectra names = {")'
printf, out_unit, format='(3(" ",a,:,","),$)', spectra_names
printf, out_unit, format='("}")'
endif
if n_elements(wave_len) ne 0 then begin
printf, out_unit, format='("wavelength = {")'
printf, out_unit, format='(7(" ",f0,:,","),$)', wave_len
printf, out_unit, format='("}")'
endif
if n_elements(wave_fwhm) ne 0 then begin
printf, out_unit, format='("fwhm = {")'
printf, out_unit, format='(7(" ",f0,:,","),$)', wave_fwhm
printf, out_unit, format='("}")'
endif
if n_elements(bad_bands) ne 0 then begin
printf, out_unit, format='("bbl = {")'
printf, out_unit, format='(10(" ",i0,:,","),$)', bad_bands
printf, out_unit, format='("}")'
endif
if n_elements(history) ne 0 then begin
printf, out_unit, format='(a)',''
printf, out_unit, format='("history = begins")'
printf, out_unit, format='(a)', history
printf, out_unit, format='("history = ends")'
endif
end
运行此代码只需要稍微熟悉 IDL 和 IDL Workbench。
2、运行准备
首先解压缩压缩文件夹(例如,H2012095004112.L1B_ISS.bz2)。 如果其他软件不容易获得,一个不错的选择是从 http://www.7-zip.org/ 免费下载 7-zip。
使用 *.h5 扩展名(例如,H2012095004112.L1B_ISS.h5)重命名生成的 HDF5 文件。 这允许 IDL 应用程序中的 HDF5 工具识别适当的格式。
如果您从本文下载了 IDL 文件,请将它们从 *.txt 重命名为 *.pro(例如,nrl_hico_h5_to_flat 到 nrl_hico_h5_to_flat.pro); 否则,如果您从 NASA 网站下载它们,它们已经具有正确的命名约定。
在 IDL Workbench 中打开 IDL 文件。 为此,只需双击文件管理器中的文件,如果您的机器上安装了 IDL,这些文件将自动在 IDL 中打开。 或者,您可以启动 ENVI+IDL 或仅启动 IDL,然后在 IDL 工作台中选择文件 > 打开。
按以下顺序编译每个文件:(i) nrl_hico_h5_to_flat.pro、(ii) byte_ordering.pro 和 (iii) write_nrl_header.pro。 在 IDL Workbench 中,这可以通过单击与给定文件关联的选项卡然后选择菜单栏中的编译按钮来实现。
您最终只会运行 nrl_hico_h5_to_flat.pro 的代码,但此应用程序依赖于其他文件; 因此,它们也需要被编译的原因。
运行 nrl_hico_h5_to_flat.pro 的代码,这是通过单击此文件的选项卡然后选择菜单栏中的运行按钮来完成的。
然后系统会提示您输入 *.h5 输入文件(例如 H2012095004112.L1B_ISS.h5),以及您希望写入输出文件的目录。
没有与此操作关联的状态栏; 但是,如果您仔细查看工作台底部 IDL 控制台中的 IDL 提示,您会注意到它会在进程运行时改变颜色,并在进程完成时恢复正常颜色。 无论如何,该过程相对较快,通常不到一分钟即可完成。
完成后,将创建两组输出文件(数据文件 + 相关的头文件),一组用于 L1B 辐射数据,一组用于导航数据。
3、HICO数据后处理
以下是准备 HICO 数据以供进一步处理所需的最后步骤:
在 ENVI 中打开 L1B 辐射和相关的导航数据。 您会注意到图像的一侧显示出包含零值的黑色条纹。
正如 HICO 网站上所述:“在 HICO 的运输和安装过程中的某个时刻,传感器相对于观察狭缝发生了物理移动。 在每个场景中都可以看到观察缝的边缘。” 通过简单地裁剪掉每个数据文件中受影响的像素,可以消除这种影响。 对于标准正向 (+XVV) 的场景,裁剪包括场景左侧的 10 个像素和右侧的 2 个像素。 相反,对于反向方向 (-XVV) 的场景,右侧裁剪 10 个像素,左侧裁剪 2 个像素。
如果您不确定特定场景的方向,可以在 hico_orientation_from_quaternion 下新创建的头文件中指定方向。
空间裁剪可以通过在 ENVI 工具箱中选择栅格管理 > 调整数据大小,选择相关的输入文件,选择空间子集的选项,将样本 11-510 的图像子集用于正向(3-502 用于反向), 并分配一个新的输出文件名。 根据需要对每个数据集重复。
HDF5 格式的 HICO 场景还需要光谱裁剪,以将总波长从 128 减少到 87 波段子集,从 0.4-0.9 um(400-900 nm)。 该子集之外的波段被认为不太准确,通常不包括在分析中。
也可以通过在 ENVI 工具箱中选择栅格管理 > 调整数据大小来执行光谱裁剪,仅在这种情况下使用光谱子集选项并选择波段 10-96(对应于 0.40408-0.89669 um)同时排除波段 1-9 和 97 -128。 此步骤只需要应用于高光谱 L1B 辐射数据。
如果需要,可以在同一步骤中应用光谱和空间裁剪。
HICO 场景现在已准备好在 ENVI 中进行进一步处理和分析。