ENVI/IDL实现HJ卫星气溶反演：http://blog.sina.com.cn/s/blog_764b1e9d01019hdw.html
这里将环境卫星气溶胶反演的三个工具和查找表建立源程序分享给大家，不同于之前的modis气溶胶反演程序，该程序做了查找表插值，因此气溶胶反演结果值是连续的。
源程序下载链接：https://pan.baidu.com/s/1oBzMI2gy_-4Cww7kyWXHIw
提取码：envi
;********************下面是角度插值工具的源码示例************************************************
;+
; :DESCRIPTION:
;
; Read HJ satelite angle data (.txt) and interpolation to polygon
;
; :AUTHOR: tiands@esrichina.com.cn
;
; :Date: 2013-7-15
;-
PRO HJ_ANGLE_EVENT,ev
COMPILE_OPT idl2
WIDGET_CONTROL, ev.TOP, get_uvalue = pstate
tagName = TAG_NAMES(ev, /structure_name)
;关闭程序事件
IF tagName EQ 'WIDGET_KILL_REQUEST' THEN BEGIN
ret = DIALOG_MESSAGE('Really Exit？', title = 'HJ Angle Tool', /question)
IF ret EQ 'No' THEN RETURN
WIDGET_CONTROL, ev.TOP, /destroy
RETURN
ENDIF
uname=WIDGET_INFO(ev.id,/uname)
CASE uname OF
;选择ccd目录
'bt_ccddir' : BEGIN
CCDPATH=DIALOG_PICKFILE(Title='Select CCD Directory:',/directory,path=pstate.Path,get_path=path)
WIDGET_CONTROL,pstate.text_ccd,set_value=CCDPATH
pstate.path=path
WIDGET_CONTROL,ev.top,set_uvalue=pstate
END
;裁剪
'evf_is' : BEGIN
WIDGET_CONTROL,pstate.BT_EVF ,sensitive=1
WIDGET_CONTROL,pstate.TEXT_EVF ,sensitive=1
pstate.isevf=1
WIDGET_CONTROL, ev.TOP, set_uvalue = pstate
END
;不裁剪
'evf_no' : BEGIN
WIDGET_CONTROL,pstate.BT_EVF ,sensitive=0
WIDGET_CONTROL,pstate.TEXT_EVF ,sensitive=0
pstate.isevf=0
WIDGET_CONTROL, ev.TOP, set_uvalue = pstate
END
;选择裁剪文件
'bt_evf' : BEGIN
evfPATH=DIALOG_PICKFILE(Title='Select Evf File:',path=pstate.Path,get_path=path,filter='*.shp')
WIDGET_CONTROL,pstate.text_evf,set_value=evfPATH
pstate.path=path
WIDGET_CONTROL,ev.top,set_uvalue=pstate
END
;选择结果保存目录
'bt_resultdir' : BEGIN
resultPATH=DIALOG_PICKFILE(Title='Select Result Directory:',/directory,path=pstate.Path,get_path=path)
WIDGET_CONTROL,pstate.text_result,set_value=resultPATH
pstate.path=path
WIDGET_CONTROL,ev.top,set_uvalue=pstate
END
;点击OK
'ok' : BEGIN
;是否裁剪判断
IF pstate.isevf EQ 1 THEN BEGIN
WIDGET_CONTROL,pstate.TEXT_EVF,get_value=shpName
shpName=shpName[0]
;shp转为evf
;prj文件路径
prjpath=FILE_DIRNAME(shpName)+'\'+STRMID(FILE_BASENAME(shpName),0,STRLEN(FILE_BASENAME(shpName))-4)+'.prj'
IF FILE_TEST(prjpath) EQ 0 THEN BEGIN
info=DIALOG_MESSAGE('不存在Prj文件！',title='提示',/cancel)
RETURN
ENDIF
evfname='C:\EvfFile.evf'
SHP2EVF,shpName,prjpath,evfname
ENDIF ELSE BEGIN
evfname=''
ENDELSE
;获取CCD目录和结果保存目录
WIDGET_CONTROL,pstate.text_ccd,get_value=ccddir
WIDGET_CONTROL,pstate.text_result,get_value=resultdir
IF (FILE_TEST(ccddir,/directory) EQ 0 ) THEN BEGIN
info=DIALOG_MESSAGE('请选择CCD文件目录！',title='提示',/cancel)
RETURN
ENDIF
IF ( FILE_TEST(resultdir,/directory) EQ 0 ) THEN BEGIN
info=DIALOG_MESSAGE('请选择结果保存目录！',title='提示',/cancel)
RETURN
ENDIF
;搜索ccd目录下的相应文件
file_tif = FILE_SEARCH(ccddir,'*.tif')
file_xml = FILE_SEARCH(ccddir,'*.xml')
file_txt = FILE_SEARCH(ccddir,'*.txt')
IF (FILE_TEST(file_tif[0]) EQ 0)THEN tif=DIALOG_MESSAGE('CCD目录下缺少TIF文件！',title='提示',/cancel)
IF (FILE_TEST(file_xml[0]) EQ 0)THEN xml=DIALOG_MESSAGE('CCD目录下缺少XML文件！',title='提示',/cancel)
IF (FILE_TEST(file_txt[0]) EQ 0)THEN txt=DIALOG_MESSAGE('CCD目录下缺少TXT文件！',title='提示',/cancel)
;;读取成像日期和时间
GETDATETIMEFROMXML, file_xml[0], otYear, otMonth, otDay, otHour, otMinute
;计算日角，赤纬角，时差
ANGLE_CAL,otYear,otMonth,otDay,sunlat,Et
;;获取文件信息
OK = QUERY_TIFF(file_tif[0], Info)
Dims = Info.DIMENSIONS
Ns = Dims[0]
Nl = Dims[1]
;临时文件路径
out_szname=resultdir+'SZ_resize.img'
out_saname=resultdir+'SA_resize.img'
out_gzname=resultdir+'GZ_resize.img'
out_ganame=resultdir+'GA_resize.img'
;这里将四个角度分开计算是因为数组占用内存太大，分开计算减小内存消耗，不至于卡机
;计算太阳天顶角
SZ_fid=SOLAR_SZ(file_tif[0],file_txt[0],Ns,Nl,evfName,out_szname,otHour,otMinute,sunlat,Et)
;计算太阳方位角
SA_fid=SOLAR_SA(file_tif[0],file_txt[0],Ns,Nl,evfName,out_saname,otHour,otMinute,sunlat,Et)
;;观测天顶角
GZ_fid=SATELLITE_GZ(file_tif[0],file_txt[0],Ns,Nl,evfName,out_gzname)
;;观测方位角
GA_fid=SATELLITE_GA(file_tif[0],file_txt[0],Ns,Nl,evfName,out_ganame)
;下面就是讲四个角度文件合成
IF (SZ_fid EQ -1 || SA_fid EQ -1 || GZ_fid EQ -1 || GA_fid EQ -1 ) THEN BEGIN
RETURN
ENDIF
ENVI_FILE_QUERY, SZ_fid,dims=dims
fid = [SZ_fid,SA_fid,GZ_fid,GA_fid]
pos = [0,0,0,0]
dims_all = LONARR(5,4) ;
FOR i=0,3 DO BEGIN
dims_all[*,i]=dims
ENDFOR
out_proj = ENVI_GET_PROJECTION(fid=SZ_fid,pixel_size=out_ps)
out_dt = 4
;将四个角度文件合成
out_name = resultdir+'Angle.dat'
out_bname=['Solar Zenith Angle','Solar Azimuth Angle','Satellite Zenith Angle','Satellite Azimuth Angle']
ENVI_DOIT, 'envi_layer_stacking_doit', $
fid=fid, pos=pos, dims=dims_all, $
out_dt=out_dt, out_name=out_name, $
interp=2, out_ps=out_ps,out_bname=out_bname, $
out_proj=out_proj, r_fid=r_fid
ENVI_FILE_MNG, id =SZ_fid,/remove,/delete
ENVI_FILE_MNG, id =SA_fid,/remove,/delete
ENVI_FILE_MNG, id =GZ_fid,/remove,/delete
ENVI_FILE_MNG, id =GA_fid,/remove,/delete
; ;关闭窗体
; WIDGET_CONTROL,ev.top,/destroy
END
'cancel' : BEGIN
WIDGET_CONTROL,ev.top,/destroy
RETURN
END
ENDCASE
END
;根据读取的角度数据，进行分块插值
PRO ANGLE_INTERP,HJCCDName,Image_angle,evfName,strHJSave,r_fid,strarr
COMPILE_OPT idl2
;将数组写出临时文件，然后获取fid和map_info
tiffile = envi_get_tmp()
ENVI_OPEN_FILE,HJCCDName,r_fid=tif_fid
map_info=ENVI_GET_MAP_INFO(fid=tif_fid)
ENVI_WRITE_ENVI_FILE,Image_angle,map_info=map_info,out_name=tiffile
ENVI_FILE_MNG,id=tif_fid,/remove
Image_angle=!null
ENVI_OPEN_FILE,tiffile,r_fid=Angle_fid
IF (FILE_TEST(evfName) EQ 1) THEN BEGIN
;在插值之前可以进行裁剪
SPATIALSUBSET,Angle_fid,evfName,resize_fid
ENVI_FILE_MNG, id =angle_fid,/remove,/delete
;这里如果返回的名字还用SZ_fid是有问题的，所以命名为resize_fid，然后再将resize_fid赋予SZ_fid
Angle_fid=resize_fid
ENDIF
;获得fid的基本信息，为分块准备
ENVI_FILE_QUERY,Angle_fid,ns=ns,sensor_type = sensor_type,$
file_type = file_type,$
nl=nl,nb=nb,dims=dims, offset = offset, bnames = bnames
pos=LINDGEN(nb)
;获取map_info用于输出
map_info=ENVI_GET_MAP_INFO(fid=Angle_fid)
;进度条
ENVI_REPORT_INIT, strarr, $
title="Get HJ Angle File!", $
base=base ,/INTERRUPT
;打开文件，并获取文件设备号
OPENW,unit,strHJSave,/get_lun
;envi分块初始化操作
tile_id=ENVI_INIT_TILE(Angle_fid,pos,num_tiles=num_tiles,xs=dims[1],$
xe=dims[2],ys=dims[3],ye=dims[4],interleave=0)
;进度条
ENVI_REPORT_INC, base, num_tiles-1
;依次读取块对应数据
FOR i=0L,num_tiles-1 DO BEGIN
;进度条
ENVI_REPORT_STAT,base, i, num_tiles-1, CANCEL=cancelvar
; 判断是否点击取消
IF cancelVar EQ 1 THEN BEGIN
tmp = DIALOG_MESSAGE('Click Cancel'+STRING(i)+'%',/info)
ENVI_REPORT_INIT, base=base, /finish
BREAK
ENDIF
;获取块数据
data=ENVI_GET_TILE(tile_id,i,band_index=band_index)
;角度插值,经过插值，这里的data变成了double类型
;所以在写出头文件的时候要注意它的data_type是5
GRIDTPS,data
;将读出的块数据存储到文件中
WRITEU,unit,data
ENDFOR
FREE_LUN,unit
ENVI_SETUP_HEAD, fname=strHJSave, $
ns=ns, nl=nl, nb=nb, $
interleave=0, data_type=4, $
offset=offset, map_info=map_info, $
bnames = bnames, $
sensor_type = sensor_type, $
file_type = file_type, $
/write, /open, r_fid = r_fid
;关闭块
ENVI_TILE_DONE,tile_id
;进度条
ENVI_REPORT_INIT, base=base, /finish
;关闭文件
ENVI_FILE_MNG,id=Angle_fid,/remove,/delete
END
;太阳天顶角
FUNCTION SOLAR_SZ,HJCCDName,angle_txt,Ns, Nl,evfName,strHJSave,otHour,otMinute,sunlat,Et
COMPILE_OPT idl2
;;创建太阳天顶角
Image_SZ = FLTARR(Ns, Nl)
;遍历角度文件
OPENR, lun, angle_txt, /get_lun
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
;;得到所有行数
iLines = FILE_LINES(angle_txt)
; iL1NlValue = LONARR(iLines - 3)
FOR ii = 0UL, iLines - 4 DO BEGIN
strTempTxt = '' & READF, lun, strTempTxt
strSplitResults = STRSPLIT(strTempTxt, /EXTRACT)
line=ROUND(FLOAT(strSplitResults[2]))-1
sample=ROUND(FLOAT(strSplitResults[3]))-1
;经纬度获取
longitude=strSplitResults[4]
Latitude=strSplitResults[5]
SOLAR_POSITION,sunlat,Et,otHour,otMinute,Latitude,longitude, otSunAA, otSunZA
;太阳天顶角ENVI有函数可以计算envi_compute_sun_angle
Image_SZ[sample,line]=otSunZA
ENDFOR
FREE_LUN, lun
;进度条STRING(13b)
strarr=['Input:'+angle_txt,'Output:'+strHJSave,$
'Processing:Solar Zenith Angle Block interpolation']
ANGLE_INTERP,HJCCDName,Image_SZ,evfName,strHJSave,r_fid,strarr
RETURN,r_fid
END
;太阳方位角计算
FUNCTION SOLAR_SA,HJCCDName,angle_txt,Ns, Nl,evfName,strHJSave,otHour,otMinute,sunlat,Et
COMPILE_OPT idl2
;;创建太阳方位角
Image_SA = FLTARR(Ns, Nl) ;
;遍历角度文件
OPENR, lun, angle_txt, /get_lun
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
;;得到所有行数
iLines = FILE_LINES(angle_txt)
; iL1NlValue = LONARR(iLines - 3)
FOR ii = 0UL, iLines - 4 DO BEGIN
strTempTxt = '' & READF, lun, strTempTxt
strSplitResults = STRSPLIT(strTempTxt, /EXTRACT)
line=ROUND(FLOAT(strSplitResults[2]))-1
sample=ROUND(FLOAT(strSplitResults[3]))-1
;经纬度获取
longitude=strSplitResults[4]
Latitude=strSplitResults[5]
SOLAR_POSITION,sunlat,Et,otHour,otMinute,Latitude,longitude, otSunAA, otSunZA
;太阳方位角
Image_SA[sample,line]=otSunAA
ENDFOR
FREE_LUN, lun
;进度条STRING(13b)
strarr=['Input:'+angle_txt,'Output:'+strHJSave,$
'Processing:Solar Azimuth Angle Block interpolation']
ANGLE_INTERP,HJCCDName,Image_SA,evfName,strHJSave,r_fid,strarr
RETURN,r_fid
END
;观测天顶角计算
FUNCTION SATELLITE_GZ ,HJCCDName,angle_txt,Ns, Nl,evfName,strHJSave
COMPILE_OPT idl2
; ;创建观测天顶角
Image_GZ = FLTARR(Ns, Nl) ;
;遍历角度文件
OPENR, lun, angle_txt, /get_lun
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
;;得到所有行数
iLines = FILE_LINES(angle_txt)
; iL1NlValue = LONARR(iLines - 3)
FOR ii = 0UL, iLines - 4 DO BEGIN
strTempTxt = '' & READF, lun, strTempTxt
strSplitResults = STRSPLIT(strTempTxt, /EXTRACT)
line=ROUND(FLOAT(strSplitResults[2]))-1
sample=ROUND(FLOAT(strSplitResults[3]))-1
; ;观测天顶角
Image_GZ[sample,line]=strSplitResults[6]
ENDFOR
FREE_LUN, lun
;进度条STRING(13b)
strarr=['Input:'+angle_txt,'Output:'+strHJSave,$
'Processing:Satellite Zenith Angle Block interpolation']
ANGLE_INTERP,HJCCDName,Image_GZ,evfName,strHJSave,r_fid,strarr
RETURN,r_fid
END
;观测方位角计算
FUNCTION SATELLITE_GA ,HJCCDName,angle_txt,Ns, Nl,evfName,strHJSave
COMPILE_OPT idl2
; ;创建观测方位角
Image_GA = FLTARR(Ns, Nl)
;遍历角度文件
OPENR, lun, angle_txt, /get_lun
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
strTempTxt = '' & READF, lun, strTempTxt
;;得到所有行数
iLines = FILE_LINES(angle_txt)
; iL1NlValue = LONARR(iLines - 3)
FOR ii = 0UL, iLines - 4 DO BEGIN
strTempTxt = '' & READF, lun, strTempTxt
strSplitResults = STRSPLIT(strTempTxt, /EXTRACT)
line=ROUND(FLOAT(strSplitResults[2]))-1
sample=ROUND(FLOAT(strSplitResults[3]))-1
; ;观测方位角
Image_GA[sample,line]=strSplitResults[7]
ENDFOR
FREE_LUN, lun
;进度条STRING(13b)
strarr=['Input:'+angle_txt,'Output:'+strHJSave,$
'Processing:Satellite Azimuth Angle Block interpolation']
ANGLE_INTERP,HJCCDName,Image_GA,evfName,strHJSave,r_fid,strarr
RETURN,r_fid
END
;;把二维离散数据插值成二维平面
;;inBlockImage：输入的离散数据，同时也是输出的二维插值结果数据
PRO GRIDTPS, inBlockImage
COMPILE_OPT idl2
;;得到插值数据的大小
Dims = SIZE(inBlockImage, /DIMENSION)
iWidth = Dims[0]
iHeight = Dims[1]
;;iIndex = Where(inBlockImage le inMax and inBlockImage ge inMin and inBlockImage ne 0, iCount)
iIndex = WHERE(inBlockImage NE 0, iCount)
;;GRID_TPS requires at least 7 noncolinear points
IF iCount LT 8 THEN BEGIN
RETURN
ENDIF
;;将一维下标转换成二维下标
NlIndex = iIndex / iWidth
NsIndex = iIndex - NlIndex * iWidth
;;插值，返回
iGoodDatas = inBlockImage[iIndex]
inBlockImage = GRID_TPS(NsIndex, NlIndex, iGoodDatas, NGRID = [iWidth, iHeight], START = [0, 0], DELTA = [1, 1])
inBlockImage=FLOAT(inBlockImage)
END
;shp转换成evf文件
PRO SHP2EVF ,shpPath,prjPath,evfPath
COMPILE_OPT IDL2
projstr=prjPath
OPENR,lun,projstr,/GET_LUN
shpPrjstr=''
READF,lun,shpPrjstr
FREE_LUN,lun
shapeproject=ENVI_PROJ_CREATE(type=42,pe_coord_sys_str=shpPrjstr)
;print,shapeproject
shapefile=shpPath
oshp=OBJ_NEW('IDLffShape',shapefile)
oshp->GETPROPERTY,n_entities=n_ent,Attribute_info=attr_info,$
n_attributes=n_attr,Entity_type=ent_type
evfpath= evfPath
evf_ptr = ENVI_EVF_DEFINE_INIT(evfpath,$
projection=shapeproject,data_type=3,$
layer_name='EVF File')
FOR i=0,n_ent-1 DO BEGIN
ent=oshp->GETENTITY(i)
;创建evf文件;当有多个polygon是，用循环
N_VERTICES=ent.N_VERTICES
parts=*(ent.PARTS);加入代码
shpsize=SIZE(parts,/N_ELEMENTS)
IF (shpsize EQ 1) THEN BEGIN
parts_ptr=[0,(N_VERTICES)];加入代码
ENDIF ELSE BEGIN
parts_ptr=[0,parts[1],(N_VERTICES)];加入代码
ENDELSE
vert=*(ent.VERTICES)
ENVI_EVF_DEFINE_ADD_RECORD,evf_ptr,vert,parts_ptr=parts_ptr,type=5
ENDFOR
evf_id = ENVI_EVF_DEFINE_CLOSE(evf_ptr, /return_id)
ENVI_EVF_CLOSE,evf_id
;attr=oshp->GetAttributes(/All)
;envi_write_dbf_file, dbfPath, attr
OBJ_DESTROY,oshp
END
;基于EVF文件的掩膜
PRO SPATIALSUBSET,angle_fid,evfName,r_fid
COMPILE_OPT IDL2
ENVI_FILE_QUERY,angle_fid,BNAMES= BNAMES,ns=ns,nl=nl,nb=nb
;打开矢量文件
evf_id = ENVI_EVF_OPEN(evfname)
;获取相关信息
ENVI_EVF_INFO, evf_id, num_recs=num_recs, $
data_type=data_type, projection=projection, $
layer_name=layer_name
roi_ids = LONARR(num_recs)
;读取各个记录的点数
FOR i=0,num_recs-1 DO BEGIN
record = ENVI_EVF_READ_RECORD(evf_id, i)
;转换为文件坐标
ENVI_CONVERT_FILE_COORDINATES,angle_fid,xmap,ymap,record[0,*],record[1,*]
;创建ROI
roi_id = ENVI_CREATE_ROI(color=4,ns = ns , nl = nl)
ENVI_DEFINE_ROI, roi_id, /polygon, xpts=REFORM(xMap), ypts=REFORM(yMap)
roi_ids[i] = roi_id
;记录XY的区间，裁剪用
IF i EQ 0 THEN BEGIN
xmin = ROUND(MIN(xMap,max = xMax))
yMin = ROUND(MIN(yMap,max = yMax))
ENDIF ELSE BEGIN
xmin = xMin <</span> ROUND(MIN(xMap))
xMax = xMax > ROUND(MAX(xMap))
yMin = yMin <</span> ROUND(MIN(yMap))
yMax = yMax > ROUND(MAX(yMap))
ENDELSE
ENDFOR
xMin = xMin >0
xmax = xMax < ns-1
yMin = yMin >0
ymax = yMax < nl-1
;创建掩膜，裁剪后掩
ENVI_MASK_DOIT,$
AND_OR =1, $
/IN_MEMORY, $
ROI_IDS= roi_ids, $ ;ROI的ID
ns = ns, nl = nl, $
/inside, $ ;区域内或外
r_fid = m_fid
out_dims = [-1,xMin,xMax,yMin,yMax]
newFile = envi_get_tmp()
ENVI_MASK_APPLY_DOIT, FID = angle_fid, POS = INDGEN(nb), DIMS = out_dims, $
M_FID = m_fid, M_POS = [0], VALUE = 0, $
OUT_BNAME= BNAMES+' Resized',IN_MEMORY=0,out_name=newFile,r_fid=r_fid
;掩膜文件ID移除
ENVI_FILE_MNG, id =m_fid,/remove
;移除evffid
ENVI_EVF_CLOSE,evf_id
;ENVI_FILE_MNG, id =angle_fid,/remove,/delete
END
;;从xml文件里面解析年月日时分秒
PRO GETDATETIMEFROMXML, inXML, otYear, otMonth, otDay, otHour, otMinute
COMPILE_OPT idl2
ON_ERROR, 2
IF FILE_TEST(inXML) EQ 0 THEN BEGIN
RETURN
ENDIF
;;打开XML文件，加载到DOM
oDocument = OBJ_NEW('IDLffXMLDOMDocument')
oDocument->LOAD, FILENAME = inXML
;;得到成像日期
oNodeList = oDocument->GETELEMENTSBYTAGNAME('imagingStartTime')
oProduct = oNodeList->ITEM(0)
oText = oProduct->GETFIRSTCHILD()
strDT = oText->GETNODEVALUE();;2009-08-01 04:41:27.13
;;分割出年、月、日、时、分、秒
iR = STRSPLIT(strDT, /EXTRACT)
iDate = STRSPLIT(iR[0], '-', /EXTRACT)
var = STRCMP(iDate,iR[0])
IF TOTAL(var) EQ 1 THEN iDate = STRSPLIT(iR[0], '/', /EXTRACT)
iTime = STRSPLIT(iR[1], ':', /EXTRACT)
otYear = iDate[0]
otMonth = iDate[1]
otDay = iDate[2]
otHour = iTime[0]
otMinute = iTime[1]
otSecond = iTime[2]
otYear = LONG(otYear)
otMonth = LONG(otMonth)
otDay = LONG(otDay)
otHour = LONG(otHour)
otMinute = LONG(otMinute)
otSecond = LONG(otSecond)
OBJ_DESTROY, oDocument
END
PRO ANGLE_CAL,otYear,otMonth,otDay,sunlat,Et
COMPILE_OPT idl2
;;日序的计算，利用年月日
month = [31,28,31 ,30,31, 30,31 ,31,30, 31,30 ,31]
dayall=0
IF(((otYear MOD 4)EQ 0) AND (( otYear MOD 100) NE 0)) OR ((otYear MOD 400) EQ 0)THEN BEGIN
run_year=29
ENDIF ELSE BEGIN
run_year=28
ENDELSE
dayall=TOTAL(month[0:otMonth-2])+otDay
;日角的计算
No=79.6764+0.2422*(otYear-1985)-FIX((otYear-1985)/4.)
;日角（弧度值）
day_angle=(2.*!PI*(dayall-No))/365.2422
;太阳赤纬角（角度值）
sunlat=0.3723+23.2567*SIN(day_angle)+0.1149*SIN(2.*day_angle)-0.1712*SIN(3.*day_angle)-0.758*COS(day_angle)+0.3656*COS(2.*day_angle)+0.0201*COS(3.*day_angle)
;时差计算
Et=0.0028-1.9857*SIN(day_angle)+9.9059*SIN(2.*day_angle)-7.0924*COS(day_angle)-0.6882*COS(2.*day_angle)
END
;;根据时间和坐标得到太阳角度
PRO SOLAR_POSITION,sunlat,Et,otHour,otMinute, Latitude, Longitude, otSunAA, otSunZA
COMPILE_OPT idl2
;北京时转换为地方时
Sd=otHour+8.+(otMinute-(120.-longitude)*4.)/60.
;时差计算
; Et=0.0028-1.9857*sin(day_angle)+9.9059*sin(2.*day_angle)-7.0924*cos(day_angle)-0.6882*cos(2.*day_angle)
;时差订正,得到真地方时
So=Sd+Et/60.
;太阳时角计算(角度值)
timeangle=(So-12.)*15.
;*太阳高度角度的计算（弧度值）
sunheightangle1 = ASIN(SIN(!pi*latitude/180.)*SIN(!pi*sunlat/180.)+COS(!pi*latitude/180.)*COS(!pi*sunlat/180.)*COS(!pi*timeangle/180.));
;弧度转角度
sunheightangle= sunheightangle1*180./!pi;
;太阳天顶角计算
sunzenith=90.-ABS(sunheightangle);
;*太阳方位角计算公式
sunposi1=ACOS(TAN(sunheightangle1)*TAN(!pi*latitude/180.)-SIN(!pi*sunlat/180.)/((COS(sunheightangle1))*(COS(!pi*latitude/180.))));
;将方位角由弧度转为角度
sunposi = sunposi1*180./!pi;
;判断是午前还是午后，因为午前午后方位角的值是不一样的,时角订正
IF(So GT 12) THEN BEGIN
sunposi = 180.+sunposi;
ENDIF ELSE BEGIN
sunposi = 180.-sunposi;
ENDELSE
otSunZA = sunzenith;
otSunAA = sunposi;
END
PRO HJ_ANGLE
COMPILE_OPT IDL2
; e=envi()
CATCH,err
IF (err NE 0) THEN BEGIN
CATCH,/cancel
!null=DIALOG_MESSAGE(!error_state.msg,/info)
RETURN
ENDIF
tlb=WIDGET_BASE(title='Get HJ-Angle File ',/column,tlb_frame_attr=1,/tlb_kill_request_events)
;选择ccd文件目录
ccd_base=WIDGET_BASE(tlb, /column, /frame,/align_center)
bt_ccddir=WIDGET_BUTTON(ccd_base,value='Select CCD File Directory',uname='bt_ccddir')
text_ccd=WIDGET_TEXT(ccd_base,xsize=50)
;设置widget_Base关键字创建单选button'
ex_base=WIDGET_BASE(tlb,/column,/frame,ysize=70)
exbase = WIDGET_BASE(ex_base,/EXCLUSIVE,/row)
eb1 = WIDGET_BUTTON(exbase,value ='裁剪',uname='evf_is')
eb2 = WIDGET_BUTTON(exbase,value ='不裁剪',uname='evf_no')
;设置默认按钮
WIDGET_CONTROL, eb2, /set_button
;evf裁剪文件
evf_base= WIDGET_BASE(ex_base,/row)
bt_evf=WIDGET_BUTTON(evf_base,value=' Shapefile',uname='bt_evf',sensitive=0)
text_evf=WIDGET_TEXT(evf_base,sensitive=0,xsize=36.5)
;结果保存目录
result_base=WIDGET_BASE(tlb, /row, /frame)
bt_resultdir=WIDGET_BUTTON(result_base,value='Output Directory',uname='bt_resultdir')
text_result=WIDGET_TEXT(result_base,xsize=30)
;确定按钮
bt_base=WIDGET_BASE(tlb,/row, /align_center)
bt_ok=WIDGET_BUTTON(bt_base,value='OK',uname='ok',xsize=60,ysize=25)
bt_cancel=WIDGET_BUTTON(bt_base,value='Cancel',uname='cancel',xsize=60,ysize=25)
;控制界面显示在屏幕中间
DEVICE, get_screen_size = screenSize
geom = WIDGET_INFO(tlb, /Geometry)
xCenter = screenSize[0] * 0.4
yCenter = screenSize[1] * 0.4
xHalfSize = geom.SCR_XSIZE / 2
yHalfSize = geom.SCR_YSIZE / 2
XOffset = (xCenter - xHalfSize)
YOffset = (yCenter - yHalfSize)
WIDGET_CONTROL, tlb,XOffset=XOffset, YOffset=YOffset
WIDGET_CONTROL, tlb,/REALIZE
pstate={text_ccd:text_ccd,$;ccd文件目录
bt_evf:bt_evf,$;选择evf文件
text_evf:text_evf,$;evf文件路径
isevf:0,$;裁剪按钮是否按下
Path:'' ,$
text_result:text_result $;结果保存目录
}
WIDGET_CONTROL,tlb,set_uvalue=pstate
;事件关联
XMANAGER,'hj_angle',tlb,/no_block
END
