【IDL程式碼庫】基於6S模型生成MODIS氣溶膠反演查詢表
阿新 • • 發佈：2022-05-27
;+ ; :Author: Hanzt ; :Email: [email protected]，歡迎討論交流 ; 更新日誌 ; 2017-04-27 新增元資訊描述 ; :Description ; 基於6S模型，構建MODIS氣溶膠反演查詢表,返回結果為8列n行的浮點型陣列 ; 第一列 太陽天頂角 ; 第二列 太陽方位角 ; 第三列 衛星天頂角 ; 第四列 衛星方位角 ; 第五列 T 大氣總透過率 ; 第六列 S 半球反照率 ; 第七列 RouA 大氣中氣溶膠反射率 ; 第八列 AOD 550nm氣溶膠光學厚度 ; ; 根據個人需求，可對輸出檔案形式進行改寫，如輸出為文字檔案，ENVI標準格式，csv檔案等 ; ; total sca :大氣整層透過率, 用T表示 ; spherical albedo :大氣半球反照率, 用S表示 ; reflectance :大氣中氣溶膠反射率,用Pa表示 ; ; 引數設定 ; exe_directory--6s.exe所在目錄，如'D:\6sWinExec' ; 研究對比發現，不同文獻、不同區域氣溶膠、太陽天頂角、方位角等資訊設定有所不同， ; ; 詳細6s引數設定參見6S使用者手冊 http://6s.ltdri.org/ ; 6s.exe及原始碼下載https://pan.baidu.com/s/1os3r0PyGgUKvowAJBBP1bA ; ; 呼叫方法--編譯後在命令列輸入以下命令： ; Build_Modis_Aerosol_Lookup_Table,ExePrograme='D:\RSdata',ThetaSs=[30],PhiRes=[60],ThetaVs=[30],PhiVs=[30] ; 其中： ; 輸入順序為 太陽天頂角，太陽方位角，衛星天頂角，衛星方位角 ; ThetaSs :太陽天頂角 ; PhiRes :太陽方位角 ; ThetaVs :衛星天頂角 ; PhiVs :衛星方位角 pro Build_Modis_Aerosol_Lookup_Table,$ AODs = AODs, $ ThetaSs = ThetaSs, $ PhiRes = PhiRes, $ ThetaVs = ThetaVs, $ Month = Month, $ Day = Day, $ IDATM = IDATM, $ IAER = IAER, $ XPS = XPS, $ IWAVE = IWAVE, $ IGROUN= IGROUN, $ output_lutFileName = output_lutFileName ;編譯規則優化 COMPILE_OPT idl2 e = envi(/current) if ~KEYWORD_SET(IDATM) then IDATM = 'No Absorption' if ~KEYWORD_SET(IAER) then IAER = 'No Aerosol' if ~KEYWORD_SET(IWAVE) then IWAVE = 'Red' if ~KEYWORD_SET(IGROUN) then IGROUN = 'VEGETA' if ~KEYWORD_SET(XPS) then XPS = 0.0 if ~KEYWORD_SET(output_lutFileName) then output_lutFileName = e.GetTemporaryFilename() ;MODIS波段代號 ;42 MODIS band 1 ( 0.6100-0.6850) ;43 MODIS band 2 ( 0.8200-0.9025) ;44 MODIS band 3 ( 0.4500-0.4825) ;45 MODIS band 4 ( 0.5400-0.5700) ;46 MODIS band 5 ( 1.2150-1.2700) ;47 MODIS band 6 ( 1.6000-1.6650) ;48 MODIS band 7 ( 2.0575-2.1825) INFO_IDATM = IDATM INFO_IAER = IAER INFO_IWAVE = IWAVE INFO_IGROUN = IGROUN INFO_XPS = XPS INFO_ThetaSs = ThetaSs INFO_ThetaVs = ThetaVs INFO_PhiRes = PhiRes INFO_AODs = AODs case IDATM of 'No Absorption':IDATM = 0 "Tropical":IDATM = 1 "Midlatitude Summer":IDATM = 2 "Midlatitude Winter":IDATM = 3 "Subartic Summer":IDATM = 4 "Subartic Winter":IDATM = 5 "US62":IDATM = 6 endcase case IAER of "No Aerosol":IAER = 0 "Continental":IAER = 1 "Maritime":IAER = 2 "Urban":IAER = 3 "Desert":IAER = 5 "Biomass":IAER = 6 "Stratospheric":IAER = 7 endcase case IWAVE of 'Red':IWAVE = 42 'Blue':IWAVE = 44 endcase case IGROUN of "VEGETA":IGROUN = 1 "CLEARW":IGROUN = 2 "SAND":IGROUN = 3 "LAKEW":IGROUN = 4 endcase XPS = 0-XPS IGEOM = 0 Month = 9 Day = 1 h = orderedHash($ 'IDATM',IDATM,$ 'IAER',IAER,$ 'IWAVE',IWAVE,$ 'IGROUN',IGROUN,$ 'XPS',XPS) print,h tic ExePrograme=filepath('6s.exe',root_dir=e.root_dir,subdirectory=['extensions','BuildModisAerosolLookupTable']) ;切換到6s.exe所在路徑 exe_directory = FILE_DIRNAME(ExePrograme) cd,exe_directory XPP = -1000 ;星測 INHOMO = 0 ;地表反射率均一地表 IDIRECT = 0 ;無方向效應 RAPP = -2 ;無大氣校正 V = 0 PhiV = 0 IGEOM = strtrim(IGEOM,1) ;衛星，自定義 MONTH = strtrim(MONTH,1) ;月份 DAY = strtrim(DAY,1) ;日期 IDATM = strtrim(IDATM,1) ;大氣模式中緯度夏季 IAER = strtrim(IAER,1) ;氣溶膠模式大陸型 V = strtrim(V,1) ;能見度 XPS = strtrim(XPS,1) ;目標物高度 XPP = strtrim(XPP,1) ;星測 IWAVE = strtrim(IWAVE,1) ;自定義1輸入波段範圍和反射相函式42為MODIS的RED INHOMO = strtrim(INHOMO,1) ;地表反射率均一地表 IDIRECT = strtrim(IDIRECT,1) ;無方向效應 IGROUN = strtrim(IGROUN,1) ;綠色植被 RAPP = strtrim(RAPP,1) ;無大氣校正 AODs = strtrim(AODs,1) ThetaSs = strtrim(ThetaSs,1) PhiRes = strtrim(PhiRes,1) ThetaVs = strtrim(ThetaVs,1) PhiV = strtrim(PhiV,1) info = strarr(1,13) in_6s = '6sIn.txt' Nlines = ulong64(AODs.length)*$ ThetaSs.length*PhiRes.length*ThetaVs.length openw,lut_lun,output_lutFileName,/GET_LUN count = 0ULL foreach AOD, AODs do $ foreach ThetaS, ThetaSs do $ foreach PhiRe, PhiRes do $ foreach ThetaV, ThetaVs do begin info[0] = IGEOM info[1] = [ThetaS+' '+PhiRe+' '+ThetaV+' '+PhiV+' '+month+' '+day] info[2] = idatm info[3] = IAER info[4] = v info[5] = AOD info[6] = xps info[7] = xpp info[8] = iwave info[9] = inhomo info[10] = idirect info[11] = igroun info[12] = rapp openw, lun_in,in_6s,/get_lun printf, lun_in, info free_lun,lun_in SPAWN,'6s.exe <6sIn.txt> 6sout.txt',/hide ;獲取6sout檔案行數，初始化一個Str用於儲存檔案內容 fileLines = file_lines('6sout.txt') Str = strarr(1,fileLines) ;開啟6sout檔案 openr,lun_out,'6sout.txt',/get_lun ;將檔案內容賦值給Str readf,lun_out,Str ;獲取T 大氣透過率 ;StartsWith僅支援IDL8.4及以後版本 !null = max(Str.StartsWith('* total sca.'),T_idx) T = strmid(Str[T_idx],17,8,/reverse_offset) ;獲取S 大氣半球反照率 !null = max(Str.StartsWith('* spherical albedo'),S_idx) S = strmid(Str[S_idx],17,8,/reverse_offset) ;獲取Pa 大氣氣溶膠反射率 !null = max(Str.StartsWith('* reflectance'),Pa_idx) RouA = strmid(Str[Pa_idx],17,8,/reverse_offset) free_lun,lun_out writeu,lut_lun,float([ThetaS,ThetaV,PhiRe,T,S,RouA,AOD]) count++ endforeach free_lun,lut_lun ; Descripe = 'Powered by Esri-China(Beijing) Modis Aerosol Lookup Table '+string(13b)+$ ; 'Using AODs of [['+AODs.join('-')+']] at 550nm'+STRING(13b)+$ ; ' IDATM of ['+STR_IDATM+']'+string(13b)+$ ; ' IAER of ['+STR_IAER+']'+string(13b)+$ ; ' IWAVE of ['+STR_IWAVE+']'+string(13b)+$ ; ' IGROUN of ['+STR_IGROUN+'] base on 6s' Descripe = 'Modis Aerosol Lookup Table Powered by Esri-China(Beijing)' ns = (float([ThetaS,ThetaV,PhiRe,T,S,RouA,AOD])).length ENVI_SETUP_HEAD,fname = output_lutFileName,$ ns = ns,nl = count,nb = 1,interleave = 0,data_type = 4,$ DESCRIP = Descripe,Bnames = 'Modis Aerosol Lookup Table',$ /write raster = e.OpenRaster(output_lutFileName) metaData = raster.metaData metaData.addItem,'Target Altitude (Km)', INFO_XPS metaData.addItem,'IDATM', INFO_IDATM metaData.addItem,'IAER', INFO_IAER metaData.addItem,'IWAVE', INFO_IWAVE metaData.addItem,'IGROUN',INFO_IGROUN metaData.addItem,'Solar Zenith', INFO_ThetaSs metaData.addItem,'Satellite Zenith', INFO_ThetaVs metaData.addItem,'Azimuthal Angle Difference', INFO_PhiRes metaData.addItem,'AODS at 550nm',INFO_AODS raster.WriteMetadata raster.close toc end