发布时间:2022-12-16 08:30
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-2allfirst_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