2. 程式人生 > 其它 >【IDL程式碼庫】環境衛星CCD資料氣溶膠反演工具原始碼分享

【IDL程式碼庫】環境衛星CCD資料氣溶膠反演工具原始碼分享

阿新 • • 發佈：2022-05-27

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: [email protected]

;

; :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

相關推薦

【IDL程式碼庫】環境衛星CCD資料氣溶膠反演工具原始碼分享

ENVI/IDL實現HJ衛星氣溶反演：http://blog.sina.com.cn/s/blog_764b1e9d01019hdw.html 這裡將環境衛星氣溶膠反演的三個工具和查詢表建立源程式分享給大家，不同於之前的modis氣溶膠反演程式，該程式做了查詢表插值，

ENVI擴充套件工具：環境衛星 CCD 資料氣溶膠反演工具

三個擴充套件工具分別實現了：計算表觀反射率、角度插值、氣溶膠反演功能。

【IDL程式碼庫】Colibri開源軟體包

Colibri介紹 Colibri是一個基於ENVI/IDL平臺實現的一個開源軟體包。裡邊包含了多種多樣的演算法和程式。

【IDL程式碼庫】IDL實現MODIS HDF檔案的輻射定標（氣溶膠系統）

一、介紹： MODIS 1KM資料的HDF包括250米和500米的兩個通道的資料，所有通道的解析度均為1000米。HDF的科學資料集包含的波段及波段型別見圖1。

【IDL程式碼庫】IDL呼叫6S.exe生成查詢表原始碼 （氣溶膠系統）

最新下載地址，包涵了原始碼和6s.exe程式：http://pan.baidu.com/s/1ppk1C ;+ ;:Description:

【IDL程式碼庫】IDL實現MODIS HDF檔案的幾何校正（氣溶膠系統）

一、介紹 MODIS HDF包含經緯度資料集，因此可以進行幾何校正。系統使用的方法是先讀取HDF中的經緯度資料集，然後建立51*51的GCP控制點（見圖1），通過對GCP點的投影轉換來對角度資料集和輻射定標後的科學資料集進行

【IDL程式碼庫】IDL實現MODIS資料雲檢測（氣溶膠系統）

一、介紹 系統採用的雲檢測演算法是多光譜綜合去雲法。該演算法從可見光反射率、紅外波段亮溫值以及亮溫差等方面綜合考慮，逐步建立一個雲檢測掩膜。對於預測的MODIS資料必須是經過了輻射定標的，熱輻射強度轉換成了

【IDL程式碼庫】利用IDLDrawWidget實現C#與IDL混合程式設計

在博文 http://blog.sina.com.cn/s/blog_764b1e9d01018v35.html 中較為詳細地介紹了IDL混合程式設計的基本手段和方法。本文以 IDLDrawWidget 為例，介紹如何進行C#與IDL之間的功能呼叫。

【IDL程式碼庫】IDL中顯示分類影象的方法

IDL作為第四代視覺化語言，具有強大的視覺化能力。根據方法的不同，可以分為四類視覺化手段，分別為快速視覺化（IDL 8.0版本開始提供）、直接圖形法、物件圖形法和智慧化程式設計工具。

【IDL程式碼庫】大資料分塊寫入HDF5檔案示例程式碼

IDL提供了專門針對HDF5科學資料格式的讀寫函式庫。可以參考IDL幫助的Routines (alphabetical) > Routines: H > HDF5 Routines 這個章節檢視函式列表。

【IDL程式碼庫】陣列編輯器

利用此程式碼可以實現手動編輯一維和二維陣列。程式碼中使用了 Clipboard 物件，支援使用者可以從其他程式（如Excel或寫字板等）複製和貼上文字。

【IDL程式碼庫】IDL中實現bwareaopen功能（過濾斑塊）

Matlab中的bwareaopen函式功能可檢視百度百科： http://baike.baidu.com/link?url=zvIOHl35oJmAcA6pSIQh1INKyfppE8IkidfVEw_p0y2JtgFz1_WKfCB1B59ZU2_3GemYozbfzNUVAL2bLVaJNq

【IDL程式碼庫】使用IDL視覺化引力波資料

新年新氣象，來點（看上去）高階的東西。 本文內容參考此文編寫（http://www.codingpy.com/article/gwpy-ligo-analyze-gravitational-waves-data/，作者：EarlGrey@程式設計派）。引力波科普知識，以及Python視覺化

【IDL程式碼庫】一個完整的ENVI擴充套件工具原始碼

以TVDI VTCI擴充套件工具為例，為廣大遙感愛好者提供一個完整ENVI/IDL二次開發示例，包括演算法編寫、資料分塊處理、繪圖、IDL介面搭建、事件響應和自定義ENVITask等內容。

【IDL程式碼庫】基於6S模型生成MODIS氣溶膠反演查詢表

;+ ; :Author: Hanzt ; :Email: [email protected]，歡迎討論交流 ; 更新日誌 ;   2017-04-27 新增元資訊描述

【IDL程式碼庫】解決IDL生成多邊形向量在ArcMap中不閉合問題

偶然發現的問題。 比如寫出一個4個頂點的矩形向量，在ArcGIS中如果設定不填充，竟然只有三個邊。。。但是在ENVI中是顯示正常的。後來請教同事後發現，是由於不閉合的原因。

【IDL程式碼庫】柵格轉向量內部介面

自ENVI 5.2版本開始提供了ENVIClassificationToShapefileTask，可以實現由分類影象轉換為shapefile向量檔案的功能。

【IDL程式碼庫】如何用IDL獲取柵格範圍平均海拔

FLAASH已提供IDL介面（http://blog.sina.com.cn/s/blog_764b1e9d0102xxrk.html），但在FLAASH大氣校正中需要輸入研究區平均海拔，我們已經介紹過如何利用ENVI的統計功能獲取研究區高程（http://blog.sina.com.cn/s/b

【IDL程式碼庫】IDL實現氣泡排序演算法

本文使用IDL實現經典的氣泡排序演算法，雖然IDL內建了排序函式sort()，但是經典的氣泡排序演算法還是值得學習。氣泡排序類似於漢諾塔遊戲，關於氣泡排序wiki百科解釋如下：

【IDL程式碼庫】我用IDL畫了個表

IDL的特長就是資料分析和視覺化，可以做很多有趣的事情。 比如本人在2014年就寫了個2048遊戲，當時確實很火，感興趣的可以去了解下。