IDL处理葵花8Himawari-8标准HSD数据——制作大气校正数据集（卫星角度数据）

IDL制作卫星天顶角和方位角数据集，对应NC文件中的SAA、SAZ

PRO temp
  COMPILE_OPT idl2
  e=envi(/headless)
  
  locRaster=e.OpenRaster('**\loc.dat')
  locData=locRaster.GetData()
  pi=3.14159265358979323846
  d2r=0.01745329252;pi/180
  r2d=57.295779513;180/pi

  cols=(locdata.dim)[0]
  lines=(locdata.dim)[1]
  outdata=make_array(cols,lines,2,/FLOAT,VALUE=0)
  lon=locData[*,*,0]
  lat=locData[*,*,1]
  alt=locData[*,*,2]#单位：km
  ;计算卫星天顶角、方位角
  satz=obslook_sat(lon, lat, alt)
  ;输出角度数据集
  outraster=ENVIRaster(outdata,URI='*\satz.dat')
  outraster.Save
  locRaster.CLose
  outraster.Close
   
  e.Close
END

代码中的’**\loc.dat’是影像的像元对应的经纬度数据，可以查看https://blog.csdn.net/qq_33339770/article/details/102957857
其中有生成loc文件的过程。

角度数据集生存代码参看：https://download.csdn.net/download/urbancorbie/10343539?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522158769212919726869052219%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.57649%2522%257D&request_id=158769212919726869052219&biz_id=1&utm_source=distribute.pc_search_result.none-task-download-2_allfirst_rank_v2~rank_v25-1
生成卫星天顶角的函数如下：

;卫星天顶角、方位角
Function obslook_sat,lon,lat,alt
  COMPILE_OPT idl2

  pi=3.14159265358979323846
  d2r=0.01745329252;pi/180
  r2d=57.295779513;180/pi
 
  cols=(lon.dim)[0]
  rows=(lon.dim)[1]
  
  clat=atan(0.993305616*tan(lat*d2r));'Rpol^2 / Req^2'=0.993305616
  clon=make_array(cols,rows,2,/FLOAT,value=d2r*lon);clon=d2r*lon
  slon=make_array(cols,rows,2,/FLOAT,value=d2r*140.7);'sub_lon'=140.7,slon=d2r*140.7
  ;计算卫星与地面距离
  ;  'Earth's polar radius'=6356.7523,'(Req^2 / Rpol^2) / Req^2'=0.00669438444
  rp=6356.7523/(sqrt(1-0.00669438444*cos(clat)*cos(clat)))+alt
  ;计算地面坐标
  x0=rp*cos(clat)*cos(clon)
  y0=rp*cos(clat)*sin(clon)
  z0=rp*sin(clat)
  ;计算卫星坐标
  ;'Distance from Earth's center to virtual satellite'=42164km
  xs=42164*cos(slon)
  ys=42164*sin(slon)
  zs=make_array(cols,rows,2,/FLOAT,value=0);zs=0
  ;计算差向量
  rx=xs-x0
  ry=ys-y0
  rz=zs-z0
  ;计算地心坐标系下的极坐标形式的差向量
  ;南方,极坐标下：x=r*sin(el)*cos(az)
  ls=sin(lat*d2r)*cos(clon)*rx+sin(lat*d2r)*sin(clon)*ry-cos(lat*d2r)*rz
  ;西方,极坐标下：y=r*sin(el)*sin(az),计算的是东方
  lse=-sin(clon)*rx+cos(clon)*ry
  ;天顶,极坐标下：z=r*cos(el)
  lz=cos(lat*d2r)*cos(clon)*rx+cos(lat*d2r)*sin(clon)*ry+sin(lat*d2r)*rz
  ;计算方位角：az=atan(y/x)
  az=atan(-lse/ls)
  az[where(ls gt 0)] += pi
  az[where(az lt 0)] += pi*2
  ;计算高度角
  rg=sqrt(rx*rx+ry*ry+rz*rz)
  el=asin(lz/rg)
  el=el*r2d
  az=az*r2d
  ;方位角修正为[-180.0,180.0]
  az[where(az GT 180.0)] -= 360.0
  ;天顶角为高度角的余角
  return,[az,90.0-el]
END