1. 程式人生 > 其它 >【IDL程式碼庫】IDL實現MODIS HDF檔案的輻射定標(氣溶膠系統)

【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。其中是MODIS第i(i=29,31,32)波段的熱輻射強度。

圖 輻射定標公式

圖3 亮溫轉換公式

二、IDL原始碼

;+

;:Description:

;    Describe the procedure.

;

;:Author: [email protected]

;

;: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