IDL处理葵花8Himawari-8标准HSD数据——制作大气校正数据集（二）
IDL制作Himawari-8太阳角度数据集（太阳方位角、天顶角）
计算卫星方位角和方位角与计算太阳方位角的原理相同，根据《卫星轨道坐标系及坐标转换》（向汝建）https://max.book118.com/html/2015/0312/13226794.shtm一文中的总结，计算方位的步骤为以下几步:
计算太阳角度，首先计算儒略日JD，需要从HSD数据中读入observation start time
openr,h8_lun,inputfile,/get_lun
;获取Observation start time
point_lun,h8_lun,46
imgtime=dblarr(1);R8
readu,h8_lun,imgtime
由于HSD数据中读入的observation start time为MJD,又知MJD=JD-2400000.5，故而
jd = imgtime[0]*1.0d + 2400000.5
imgtime要时一个后面带d的double型数，乘以1.0d就可以使数据转为double类型，否则IDL会出现搞笑一幕：
如果是UTC时间，在IDL,Python或者matlab中都有将UTC时间转为儒略日时间的函数，可直接转换。
JD计算完成后是计算观测点的ECI坐标，即笛卡尔坐标系下的坐标。
计算ECI坐标，这里在之前的生成Himawari-8几何校正文件一文中https://blog.csdn.net/qq_33339770/article/details/102957857，对影像上每个像元的经纬度都进行过计算，故可形成像元坐标文件，这里的坐标是lon,lat,alt（可用ENVI自带DEM数据提取高程，也可直接定义为0）。
也可以自己定义左上角坐标80.0，60.0，行列数6001或12001，像元宽高0.02或0.01来书写lon，lat文件，但是最终计算结果会在纬度为0.0的一行结果反向，不想单独处理这一行，故而选择使用通过HSD文件中的行列偏移计算出的坐标文件。
需要将lon,lat,alt转为ECI坐标系下的x0,y0,z0坐标，根据论文中所述，地理经度在ECI坐标中是时间的函数，故而首先计算地方恒星时:地方恒星时=Greenwich恒星时+lon
Greenwich恒星时的查找表比较麻烦，可以选择自己计算，计算的方法在论文中描述如下：
首先计算UTC时间为0：00的恒星时：
ut是这一日的儒略日小数，也就是上图公式中的ΔT，由于儒略日的起始时间是公元前4713年1月1日中午12点，故而回溯半天（+0.5）再对1取余。
根据论文计算UTC 0：00的恒星时为（单位为秒）：
UT = (jd+0.5) mod 1.0
jd_2000 = 2451545.0
t_cen = (jd - UT - jd_2000) / 36525.0
GMST = 24110.54841 + t_cen * (8640184.812866 + t_cen *(0.093104 - t_cen * 6.2E-6))
然后计算ΔT时间的恒星时，并将恒星时化为正值：
omega_e = 1.00273790934; Earth rotations per sidereal day
seconds_per_day = 86400.0
GMST = ( GMST + seconds_per_day * omega_E * UT) MOD seconds_per_day
GMST[where(GMST LT 0.0)] += seconds_per_day
omega_e为地球每恒星日转过的弧度
之后将恒星时转为弧度：
rval = 2.0 * pi * GMST / seconds_per_day
整段计算恒星时的代码如下(此段theta_JD函数参考https://github.com/larrylart/Unimap/blob/db5c5b503d441bb784c7ed551741233ac7350e08/src/sky/sat_code/observe.cpp中的函数ThetaG):
;计算格林威治时间并转为角度
Function thetaG_JD,jd
COMPILE_OPT idl2
pi=3.14159265358979323846
omega_e = 1.00273790934; Earth rotations per sidereal day
UT = (jd+0.5) mod 1.0
seconds_per_day = 86400.0
jd_2000 = 2451545.0
t_cen = (jd - UT - jd_2000) / 36525.0
GMST = 24110.54841 + t_cen * (8640184.812866 + t_cen *(0.093104 - t_cen * 6.2E-6))
GMST = ( GMST + seconds_per_day * omega_E * UT) MOD seconds_per_day
GMST[where(GMST LT 0.0)] += seconds_per_day
rval = 2.0 * pi * GMST / seconds_per_day
return,rval
END
将地理经度带入公式，其后即可计算地方恒星时
d2r=0.01745329252;pi/180
theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)
有了地方恒星时，则可以将lon,lat,alt转为ECI坐标下的x0,y0,z0:
;经纬度转球面坐标(WGS_72椭球)
Function ll2xyz_ellips,jd,lon,lat,alt
COMPILE_OPT idl2
pi=3.14159265358979323846
d2r=0.01745329252;pi/180
r2d=57.295779513;180/pi
Re=6378.14;赤道半径6378.140km(长半轴)
theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)
;转换到ECI坐标系下
f=1.0/298.256;地球扁率
c=1/sqrt(1+f*(f-2)*sin(lat*d2r)*sin(lat*d2r))
S=(1-f)*(1-f)*C
;计算ECI坐标系下的坐标x,y,z
x=(Re+alt/1000.0)*C*cos(lat*d2r)*cos(theta)
y=(Re+alt/1000.0)*C*cos(lat*d2r)*sin(theta)
z=(Re+alt/1000.0)*S*sin(lat*d2r)
return,hash('x0',x,'y0',y,'z0',z)
END
上述代码将经纬度坐标转为椭球体上的xyz坐标，也可以将地球视为球体，转为球面的坐标：
;经纬度转球面坐标(球体)
Function ll2xyz_sphere,jd,lon,lat,alt
COMPILE_OPT idl2
pi=3.14159265358979323846
d2r=0.01745329252;pi/180
r2d=57.295779513;180/pi
Re=6378.14;赤道半径6378.140km(长半轴)
theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)
x=Re*cos(lat*d2r)*cos(theta)
y=Re*cos(lat*d2r)*sin(theta)
z=Re*sin(lat*d2r)
return,hash('x0',x,'y0',y,'z0',z)
END
对于ECI坐标计算存疑的，可以带入论文中的一套数据来验证计算出的x0,y0,z0是否正确,作者验证过椭球ECI坐标为正确的。
有了观测点的ECI坐标（x0,y0,z0）,同时可从HSD数据集中读到Sun’s position(ECI坐标)
;获取Sun's position
point_lun,h8_lun,510
sun_pos=dblarr(3);R8
readu,h8_lun,sun_pos
则可以计算观测点与太阳的向量差[rx，ry，rz]=[xs-x0，ys-y0，zs-z0]
这个向量差为ECI坐标系下的差向量，需要转为测站坐标系（地心坐标系）下的极坐标形式的向量差，就可以得到方位角、高度角（与天顶角互余）。
这里θ为天顶角（高度角的余角），φ为方位角
计算太阳方位角、天顶角的代码如下：
;太阳天顶角、方位角
Function obslook_solar,jd,sun_pos,lon,lat,alt
COMPILE_OPT idl2
pi=3.14159265358979323846
d2r=0.01745329252;pi/180
r2d=57.295779513;180/pi
Re=6378.14;赤道半径6378.140km
;计算经纬度对应的ECI坐标
xyz=ll2xyz_ellips(jd,lon,lat,alt)
theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)
x0=xyz['x0']
y0=xyz['y0']
z0=xyz['z0']
OBJ_DESTROY,xyz
;计算相对距离
xs=sun_pos[*,*,0]
ys=sun_pos[*,*,1]
zs=sun_pos[*,*,2]
;卫星的位置是ECI坐标系，定义为[xs，ys，zs]，观察者是[x0，y0，z0]
;那么距离向量就是[rx，ry，rz]=[xs-x0，ys-y0，zs-z0]
rx=xs-x0
ry=ys-y0
rz=zs-z0
;[rx，ry，rz]是ECI坐标,要生成视角,需要转为地心地平线系统
;这个系统的z轴指向天顶,x轴指向南,y轴指向西(地心坐标),地心坐标系随地球旋转
ls=sin(lat*d2r)*cos(theta)*rx+sin(lat*d2r)*sin(theta)*ry-cos(lat*d2r)*rz;南方,极坐标下：x=r*sin(el)*cos(az)
lse=-sin(theta)*rx+cos(theta)*ry;西方,极坐标下：y=r*sin(el)*sin(az)
lz=cos(lat*d2r)*cos(theta)*rx+cos(lat*d2r)*sin(theta)*ry+sin(lat*d2r)*rz;天顶,极坐标下：z=r*cos(el)
;az:AzimuthAngle,方位角
az=atan(-lse/ls);atan(y/x)
az[where(ls gt 0)] += pi
az[where(az lt 0)] +=pi*2
;相对距离
rg=sqrt(rx*rx+ry*ry+rz*rz)
;el:ZenithAngle,天顶角
el=asin(lz/rg)
;弧度转角度
el=el*r2d
az=az*r2d
;输出SOA,SOZ
az[where(az GT 180.0)] -= 360.0
return,hash('SOA',az,'SOZ',90.0-el)
END
对角度的计算验证也可以带入论文中的一组数据，这里作者未带入数据，直接拿结果与NC文件中的角度数据集对比正确。
代码传入数据说明，此文中的imgtime）,lon,lat,alt均与矫正裁剪后影像同等大小，为60016001或1200112001大小的数组，sun_pos需要存储三维的太阳坐标（xs,ys,zs）故数组大小为600160013或12001120013
imgtime数组的数值必须带d的double型数，包装成数组时可以采取下面的方法：
obstime=make_array(cols,rows,value=imgtime[0]*1.0d)
这样的话JD的计算就不需要*1.0d
jd = obstime+ 2400000.5
由于HSD每个波段分为10个分块，每个分块的Observation start time和Sun’s position都有变化，为准确的计算太阳角度，需要将10块imgtime和sun_pos读取之后作为此波段的time和sunpos数组，之后与影像一起进行几何校正和裁剪，即可作为角度数据集的输入数组。