【IDL代码库】IDL实现MODIS HDF文件的辐射定标(气溶胶系统)

一、介绍:

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。其中【IDL代码库】IDL实现MODIS <wbr>HDF文件的辐射定标(气溶胶系统)是MODIS第i(i=29,31,32)波段的热辐射强度。

图 辐射定标公式

图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

posted @ 2022-05-27 11:47  ENVI-IDL技术殿堂  阅读(676)  评论(0编辑  收藏  举报