一、介绍:
MODIS 1KM数据的HDF包括250米和500米的两个通道的资料,所有通道的分辨率均为1000米。HDF的科学数据集包含的波段及波段类型见图1。
图1 HDF数据集
气溶胶系统中利用了影像的部分波段,因此为了节省处理的时间把波段1、2、3、6、7、8、26、29、31、32单独进行辐射定标,然后合成为一个新的数据。辐射定标公式见图2,参数可以从HDF中的数据集中读取。辐射定标之后太阳反射波段1、2、3、6、7、8、26定标后为发射率,热辐射波段29、31、32定标后辐射亮度(热辐射强度),辐射亮度转换成亮温的公式见图3。其中是MODIS第i(i=29,31,32)波段的热辐射强度。
图 2 辐射定标公式
图3 亮温转换公式
二、IDL源码
;+
;:Description:
; Describe the procedure.
;
;:Author: tiands@esrichina.com.cn
;
;:Date: 2013-9-5 11:10:27
;-
;;;=======================================================================
;;; 该功能把hdf数据中的波段1、2、3、4、6、7、8、19、26、29、31、32转换成反射率和亮温以
;;;便进行 反演
;;;======================================================================
PRO RADIATION_CORRECTION,modisname,outname
COMPILE_OPT idl2
; ENVI, /RESTORE_BASE_SAVE_FILES
; ENVI_BATCH_INIT, LOG_FILE = 'batch_log.txt'
;科学数据读取
ref250=READ_DATASET(modisname,'EV_250_Aggr1km_RefSB');读取band1 and band2
ref500=READ_DATASET(modisname,'EV_500_Aggr1km_RefSB');读取band3-band7
ref1000=READ_DATASET(modisname,'EV_1KM_RefSB') ;读取band8-band19 和band26
Emi1000=READ_DATASET(modisname,'EV_1KM_Emissive') ;读取band20-band25 和 band27-band36
;;;=======================================================================
;;; 科学数据合成这里选择了本系统中用到的波段
;;;======================================================================
ENVI_WRITE_ENVI_FILE,ref250[*,*,0],r_fid=fid_ref250_1,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref250[*,*,1],r_fid=fid_ref250_2,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref500[*,*,0],r_fid=fid_ref500_3,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref500[*,*,1],r_fid=fid_ref500_4,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref500[*,*,3],r_fid=fid_ref500_6,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref500[*,*,4],r_fid=fid_ref500_7,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref1000[*,*,0],r_fid=fid_ref1000_8,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref1000[*,*,13],r_fid=fid_ref1000_19,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,ref1000[*,*,14],r_fid=fid_ref1000_26,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,Emi1000[*,*,8],r_fid=fid_Emi1000_29,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,Emi1000[*,*,10],r_fid=fid_Emi1000_31,/IN_MEMORY
ENVI_WRITE_ENVI_FILE,Emi1000[*,*,11],r_fid=fid_Emi1000_32,/IN_MEMORY
;;;=======================================================================
;;; 传入hdf文件波段数据以及数据集索引和对应数据集属性索引(以0开始),传入的数字具体
;含义可以使用envi打开hdf的数据集,对照着数字可以找到所要的信息
;;;======================================================================
;EV_250_Aggr1km_RefSB数据集反射率
r_fid_ref250_1=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref250_1,6,8,9,0,0)
r_fid_ref250_2=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref250_2,6,8,9,1,1)
;EV_500_Aggr1km_RefSB数据集反射率
r_fid_ref500_3=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref500_3,9,8,9,0,0)
r_fid_ref500_4=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref500_4,9,8,9,1,1)
r_fid_ref500_6=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref500_6,9,8,9,3,3)
r_fid_ref500_7=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref500_7,9,8,9,4,4)
;EV_1KM_RefSB数据集校正为反射率
r_fid_ref1000_8=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref1000_8,2,8,9,0,0)
r_fid_ref1000_19=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref1000_19,2,8,9,13,13)
r_fid_ref1000_26=REFLECT_RADIANCE_CALCULATE(modisname,fid_ref1000_26,2,8,9,14,14)
;EV_1KM_Emissive数据集校正为辐射亮度
r_fid_Emi1000_29=REFLECT_RADIANCE_CALCULATE(modisname,fid_Emi1000_29,4,5,6,8,8)
r_fid_Emi1000_31=REFLECT_RADIANCE_CALCULATE(modisname,fid_Emi1000_31,4,5,6,10,10)
r_fid_Emi1000_32=REFLECT_RADIANCE_CALCULATE(modisname,fid_Emi1000_32,4,5,6,11,11)
;传入的数据是根据公式 , Ki,1=C1/r5,其中r为波长,C1为固定值,可参阅辐射亮度转换成亮温
bt_fid_29=BRIGHT_TEMPERATURE(r_fid_Emi1000_29,'1682.770175/alog(1+2606.736131/b1)')
bt_fid_31=BRIGHT_TEMPERATURE(r_fid_Emi1000_31,'1304.413871/alog(1+729.541636/b1)')
btbt_fid_32=BRIGHT_TEMPERATURE(r_fid_Emi1000_32,'1196.978785/alog(1+474.684780/b1)')
ENVI_FILE_QUERY, r_fid_ref250_1,dims=dims
;将所有科学数据合成一个图像,便于后续的几何校正工作
OUT_BNAME=['1KM Reflectance (band1) [250 Aggr]','1KM Reflectance (band2) [250 Aggr]',$
'1KM Reflectance (band3) [500 Aggr]','1KM Reflectance (band4) [500 Aggr]','1KM Reflectance (band6) [500 Aggr]',$
'1KM Reflectance (band7) [500 Aggr]','1KM Reflectance (band8)','1KM Reflectance (band19)','1KM Reflectance (band26)',$
'1KM Emissive (band29)','1KM Emissive (band31)','1KM Emissive (band32)']
ENVI_DOIT, 'cf_doit',fid=[r_fid_ref250_1,r_fid_ref250_2,r_fid_ref500_3,r_fid_ref500_4,r_fid_ref500_6,r_fid_ref500_7,$
r_fid_ref1000_8,r_fid_ref1000_19,r_fid_ref1000_26,bt_fid_29,bt_fid_31,btbt_fid_32], $
pos=[0,0,0,0,0,0,0,0,0,0,0,0], dims=dims, OUT_BNAME=OUT_BNAME,$
remove=0,r_fid=data_fid,out_name=outname
; ENVI_BATCH_EXIT
END
FUNCTION READ_DATASET,filename,dataset
SD_id = HDF_SD_START(filename)
index = HDF_SD_NAMETOINDEX(SD_id,dataset)
sds_id=HDF_SD_SELECT(SD_id,index)
HDF_SD_GETDATA,sds_id,data
HDF_SD_ENDACCESS,sds_id
RETURN,data
END
FUNCTION REFLECT_RADIANCE_CALCULATE,modisname,fid,dataset_number,scalesattr_index,offsetsattr_index,scales_index,offsets_index
COMPILE_OPT idl2
;选择数据集
SD_id = HDF_SD_START(modisname,/Read)
SdsID=HDF_SD_SELECT(SD_id,dataset_number)
;得到属性名称和属性值
HDF_SD_ATTRINFO,SdsID,scalesattr_index,NAME=re_scales,DATA=scalesData
HDF_SD_ATTRINFO,SdsID,offsetsattr_index,NAME=re_offsets,DATA=offsetsData
;反射率计算公式 R=reflectance_scale*(DN-reflectance _offset)
;辐射亮度计算公式 R=radiance_scale*(DN-radiance _offset)
ENVI_FILE_QUERY, fid, dims=dims
t_fid = [fid]
pos = [0]
;表达式
exp=STRTRIM(scalesData[scales_index],2)+'*[b1-('+STRTRIM(offsetsData[offsets_index],2)+')]'
ENVI_DOIT, 'math_doit', $
fid=t_fid, pos=pos, dims=dims,$
exp=exp,r_fid=r_fid, /IN_MEMORY
ENVI_FILE_MNG, id=fid, /remove
RETURN,r_fid
END
FUNCTION BRIGHT_TEMPERATURE,bt_fid,exp
;辐射亮度转换为亮温公式T=K2/ln(1+K1/R)其中对于特定波段K1K2是固定的
;对于band31 K1=729.541636 K2=1304.413871
; band29 K1=2606.736131 K2=1682.770175
ENVI_FILE_QUERY, bt_fid, dims=dims
t_fid = [bt_fid]
pos = [0]
;表达式
exp=exp
ENVI_DOIT, 'math_doit', $
fid=t_fid, pos=pos, dims=dims,$
exp=exp,r_fid=r_fid, /IN_MEMORY
ENVI_FILE_MNG, id=bt_fid, /remove
RETURN,r_fid
END