IDL调用6S大气校正

sixs为6s源码编译后的exe全路径

sixs='*********\sixsV2.1.exe'

代码传入的参数中有paraID，表示了影像是否为纯水域（完全在海上）

Function sixs_atmos,sixs,temp_path,Parameter
  COMPILE_OPT idl2
  ;Parameter=[SOZ,SOA,SAZ,SAA,MONTH,DAY,paraID,xps,iwave,idatm,rapp]
  SOZ=Parameter[0];太阳天顶角
  SOA=Parameter[1];太阳方位角
  SAZ=Parameter[2];卫星天顶角
  SAA=Parameter[3];卫星方位角
  MONTH=uint(Parameter[4]);月
  DAY=uint(Parameter[5]);日
  paraID=uint(Parameter[6]);1-有陆地 2-纯水域
  xps=Parameter[7];高程
  print,'平均高程',xps
  iwave=uint(Parameter[8]);光谱响应函数
  idatm=uint(Parameter[9]);大气模式
  ;rapp>0：输入为大气层顶辐射亮度
  ;-1.0<rapp<0 激活，输入为表观反射率，确保表观反射率为正值
  rapp=Parameter[10]
    
  
  ;生成对应波段的int.txt文档
  inname=temp_path+'input.txt'
  openw,lun,inname,/get_lun
  igeom=0;几何条件，自定义
  printf,lun,igeom
  printf,lun,SOZ,SOA,SAZ,SAA,MONTH,DAY;太阳天顶角/方位角 卫星天顶角/方位角
  ;2.大气模式,根据经纬度确定
  printf,lun,idatm
  ;3、气溶胶模式
  if paraID eq 2 then iaer=2;海洋型
  if paraID eq 1 then iaer=1;大陆型
  printf,lun,iaer
  ;4、气溶胶含量，使用能见度
  if paraID eq 2 then v=40;海洋型
  if paraID eq 1 then v=15;大陆型
  printf,lun,v
  ;5、目标海拔高度
  printf,lun,xps
  ;6、传感器海拔高度
  xpp=-1000
  printf,lun,xpp
  ;7、光谱参数
  printf,lun,iwave
  ; (8)地表反射率类型（ground reflectance type）,参数名称：inhomo    取值范围：0，1
  ;参数名称：inhomo    取值范围：0，1
  ; 下垫面非均一、朗伯体
  inhomo=0;均匀表面
  printf,lun,inhomo
  idirec=0;无方向效应
  printf,lun,idirec
  ;有陆地内陆用湖水 4 lake water
  if paraID eq 2 then igroun=2
  ,纯水域用纯水2，clear water
  if paraID eq 1 then igroun=4
  printf,lun,igroun
  ;9、大气订正方式
  irapp=0;-1不进行大气订正，0 朗伯体和idirec对应，1 BRDF各项异质性
  printf,lun,irapp
  ;rapp=-0.5;表观反射率大气校正
  ;rapp=0.5;辐射亮度大气校正
  printf,lun,rapp
  ;极化方式
  irop=4;1-3极化方式不同，其他值不考虑极化方式
  printf,lun,irop
  FREE_LUN,lun

  ;参数写入完毕，调用6s.exe,输出txt
  outtxtname=temp_path+"out.txt"
  ;调用6s <>符号为dos系统下的重定向符号    <从文件读取命令输入>将输出结果写入文件
  SPAWN,sixs + '<'+ inname + '>' + ' ' + outtxtname,/hide
  ;输出文本中找到xa xb xc
  nLines=file_lines(outtxtname)
  openr,lun,outtxtname,/GET_LUN
  result_arr=STRARR(nLines);行数
  readf,lun,result_arr
  FREE_LUN,lun
  ;获取6s计算结果out.txtx中的参数 
  ;辐射亮度参数
  ;'* coefficients xa xb xc : 0.00272 0.11681 0.17263 *'
  if rapp gt 0 then begin
    result_index=163
    splitxt='* coefficients xa xb xc :'
  endif
  ;表观反射率参数
  ;'* coefficients xap xb xc 0.867103 -0.209570 0.172645 *'
  if rapp gt -1 and rapp lt 0 then begin
    result_index=165
    splitxt='* coefficients xap xb xc :'
  endif
  ;提取参数
  result_str=strsplit(result_arr[result_index],splitxt,/extract)
  ;结果
  result_number=float(result_str)
  return,result_number
END