+
; :Author: Hanzt
; :Email: hanzt@esrichina.com.cn，欢迎讨论交流
; 更新日志
; 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
