2. 程式人生 > 其它 >【IDL程式碼庫】IDL實現MODIS HDF檔案的幾何校正（氣溶膠系統）

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

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

一、介紹

MODIS HDF包含經緯度資料集，因此可以進行幾何校正。系統使用的方法是先讀取HDF中的經緯度資料集，然後建立51*51的GCP控制點（見圖1），通過對GCP點的投影轉換來對角度資料集和輻射定標後的科學資料集進行幾何校正。在幾何校正的過程中將角度資料重取樣到1000*1000解析度，以便後來進行氣溶膠反演的時候，按行列號查詢角度資料。

二、IDL原始碼

原始碼下載：http://pan.baidu.com/s/1dDcXOEP

;+

;:Description:

;    Describe the procedure.

;

;:Author: [email protected]

;

;:Date: 2013-9-5 11:10:27

;-

;;=======================================================================  ;;;  改程式實現對輻射校正後的科學資料幾何校正以及對角度資料校正，輸出控制點檔案  ;;;======================================================================  PRO MODIS_REGISTER,filename,fs_name,outname_angle,outname_data,GCP1,GCP2    ENVI, /RESTORE_BASE_SAVE_FILES
    ENVI_BATCH_INIT, LOG_FILE = 'batch_log.txt'    COMPILE_OPT idl2        ;  filename ='H:\MODIS\modis資料\1km\MOD021KM.A2012126.0245.005.2012128183247.hdf'    ;  fs_name='C:\Users\tiandeshan\Desktop\idl\fushejiaozheng.img'    ;  outname_angle='C:\Users\tiandeshan\Desktop\idl\angledata.img'    ;  outname_data='C:\Users\tiandeshan\Desktop\idl\scaledata.img'
    ;  GCP1='C:\Users\tiandeshan\Desktop\idl\GCP1.dat'    ;  GCP2='C:\Users\tiandeshan\Desktop\idl\GCP2.dat'            ;經緯度讀取   Lon=READ_DATASET(filename,'Longitude')   Lat=READ_DATASET(filename,'Latitude')    ;重取樣經緯度   ENVI_WRITE_ENVI_FILE,Lon[*,*],r_fid=fid_Lon,out_name='c:/temp'    FIDlon= MODIS_RESIZE(fid_Lon)    ENVI_FILE_QUERY, FIDlon,dims=dimslon    datalon = ENVI_GET_DATA(fid=FIDlon, dims=dimslon, pos=0)       ENVI_WRITE_ENVI_FILE,Lat[*,*],r_fid=fid_Lat,out_name='c:/temp'    FIDlat= MODIS_RESIZE(fid_Lat)    ENVI_FILE_QUERY, FIDlat,dims=dimslat    datalat = ENVI_GET_DATA(fid=FIDlat, dims=dimslat, pos=0)   ;;;=======================================================================    ;;;  GCP控制點建立   ;;;======================================================================   OPENW,lun,GCP1,/Get_lun    PRINTF,lun,' gcp[0,l] ','gcp[1,l] ','gcp[2,l]  ','gcp[2,l]  ','n  ','m  '    length=0L    ;定義控制點為51*51    gcp= DBLARR(4,2601)    rstr = ["Input File :" , "Output File :"]   ENVI_REPORT_INIT,rstr,title="output the GCP", base=base,/INTERRUPT    FOR i=3.5d,2030,40.53d DO BEGIN      FOR j=3.5d,1354,27.01d DO BEGIN       gcp[0,length]=datalon[j,i]       gcp[1,length]=datalat[j,i]       gcp[2,length]=j       gcp[3,length]=i       PRINTF,lun,gcp[0,length],gcp[1,length],gcp[2,length],gcp[3,length]       length=length+1L     ENDFOR     ENVI_REPORT_STAT, base, i*1354, 2601    ENDFOR    ENVI_REPORT_INIT, base=base, /finish    FREE_LUN,lun   ;;;=======================================================================    ;;;  GCP控制點投影轉換   ;;;======================================================================   OPENW,lun,GCP2,/Get_lun   PRINTF,lun,'ogcp[0,l]','ogcp[1,l] ','ogcp[2,l]','ogcp[2,l]'    iproj= ENVI_PROJ_CREATE(/geographic)    units = envi_translate_projection_units('Meters')    datum = 'WGS-84'    oproj= ENVI_PROJ_CREATE(/utm, zone=51,datum=datum,units=units)   ENVI_CONVERT_PROJECTION_COORDINATES, gcp[0,*], gcp[1,*],iproj,outx, outy, oproj    ;    print,gcp[0,*]    ogcp= DBLARR(4,2601)    length=0L    FOR i=3.5d,2030,40.53d DO BEGIN      FOR j=3.5d,1354,27.01d DO BEGIN       ogcp[0,length]=outx[[(i-3.5)/40.53]*51+[(j-3.5)/27.01]]       ogcp[1,length]=outy[[(i-3.5)/40.53]*51+[(j-3.5)/27.01]]       ogcp[2,length]=gcp[2,length]       ogcp[3,length]=gcp[3,length]       PRINTF,lun,ogcp[0,length],ogcp[1,length],ogcp[2,length],ogcp[3,length]       length=length+1L     ENDFOR    ENDFOR    FREE_LUN,lun   ;;;=======================================================================    ;;; 開始校正   ;;;======================================================================    ;角度校正   OUT_BNAME=['Warp(SensorZenith)','Warp(SensorAzimuth)',$     'Warp(SolarZenith)','Warp(SolarAzimuth)']   anglefid=ANGLE_DATA(filename)   ENVI_FILE_QUERY,anglefid,nb=nb_angle,dims=dims_angle   REGISTER,anglefid,nb_angle,dims_angle,ogcp,oproj,outname_angle,OUT_BNAME    ;科學資料校正   OUT_BNAME=['Warp(1KM Reflectance (band1) [250 Aggr])','Warp(1KM Reflectance (band2) [250 Aggr])',$     'Warp(1KM Reflectance (band3) [500 Aggr])','Warp(1KM Reflectance (band4) [500 Aggr])','Warp(1KM Reflectance (band6) [500 Aggr])',$     'Warp(1KM Reflectance (band7) [500 Aggr])','Warp(1KM Reflectance (band8))','Warp(1KM Reflectance (band19))','Warp(1KM Reflectance (band26))',$     'Warp(1KM Emissive (band29))','Warp(1KM Emissive (band31))','Warp(1KM Emissive (band32))']           ;資料讀取    IF (ENVI_IS_GDB(fs_name)) THEN BEGIN     ENVI_OPEN_GDB,fs_name,r_fid=data_fid    ENDIF ELSE BEGIN     ENVI_OPEN_FILE,fs_name,r_fid=data_fid    ENDELSE    ;envi_open_file, fs_name, r_fid=data_fid    ENVI_FILE_QUERY, data_fid, dims=dims_data, nb=nb_data   REGISTER,data_fid,nb_data,dims_data,ogcp,oproj,outname_data,OUT_BNAME    ENVI_BATCH_EXIT  END        FUNCTION  ANGLE_DATA,filename    COMPILE_OPT idl2    ;角度資訊讀取   sensorzenith=READ_DATASET(filename,'SensorZenith')   sensorazimuth=READ_DATASET(filename,'SensorAzimuth')   solarzenith=READ_DATASET(filename,'SolarZenith')   solarazimuth=READ_DATASET(filename,'SolarAzimuth')   ;;;=======================================================================    ;;;  角度資訊合成   ;;;======================================================================   ENVI_WRITE_ENVI_FILE,sensorzenith[*,*],r_fid=fid_sez,/IN_MEMORY   ENVI_WRITE_ENVI_FILE,sensorazimuth[*,*],r_fid=fid_sea,/IN_MEMORY   ENVI_WRITE_ENVI_FILE,solarzenith[*,*],r_fid=fid_soz,/IN_MEMORY   ENVI_WRITE_ENVI_FILE,solarazimuth[*,*],r_fid=fid_soa,/IN_MEMORY   ENVI_FILE_QUERY,fid_sez,dims=sezdims   fid_sez_c=ANGLE_CALCULATE(fid_sez)   fid_sea_c=ANGLE_CALCULATE(fid_sea)   fid_soz_c=ANGLE_CALCULATE(fid_soz)   fid_soa_c=ANGLE_CALCULATE(fid_soa)   ;將四個有用資訊的影象合成一個影象，便於後續的幾何校正工作    ENVI_DOIT, 'cf_doit',fid=[fid_sez_c,fid_sea_c,fid_soz_c,fid_soa_c], $     pos=[0,0,0,0], dims=sezdims, $     remove=0,r_fid=info_fid ,/in_memory    ;將四個單獨的幾何資訊檔案釋放    ENVI_FILE_MNG, id=fid_sez, /remove    ENVI_FILE_MNG, id=fid_sea, /remove    ENVI_FILE_MNG, id=fid_soz, /remove    ENVI_FILE_MNG, id=fid_soa, /remove    angle_fid= MODIS_RESIZE(info_fid)        RETURN,angle_fid  END      FUNCTION ANGLE_CALCULATE,fid    COMPILE_OPT idl2        ENVI_FILE_QUERY, fid, dims=dims    t_fid = [fid]    pos = [0]    ;表示式    exp='b1*0.0100'    ENVI_DOIT, 'math_doit', $     fid=t_fid, pos=pos, dims=dims,$     exp=exp,r_fid=r_fid, /IN_MEMORY    ENVI_FILE_MNG, id=fid, /remove    RETURN,r_fid  END          FUNCTION READ_DATASET,filename,dataset    SD_id = HDF_SD_START(filename)    index = HDF_SD_NAMETOINDEX(SD_id,dataset)   sds_id=HDF_SD_SELECT(SD_id,index)   HDF_SD_GETINFO,sds_id,DIMS=dim   HDF_SD_GETDATA,sds_id,data   HDF_SD_ENDACCESS,sds_id    RETURN,data  END        FUNCTION MODIS_RESIZE,fid    COMPILE_OPT IDL2    ; Open the input file    IF (fid EQ -1) THEN BEGIN     RETURN,-1    ENDIF    ENVI_FILE_QUERY, fid, dims=dims, nb=nb    pos = LINDGEN(nb)    out_name = 'testimg'    ;把資料集重取樣成1354*2030    ENVI_DOIT, 'resize_doit', $     fid=fid, pos=pos, dims=dims, $     interp=1, rfact=[.20015,.2], $     out_name=out_name, r_fid=r_fid    RETURN,r_fid  END    PRO REGISTER,fid,nb,dims,ogcp,oproj,out_name,OUT_BNAME    COMPILE_OPT idl2    pos = LINDGEN(nb)    pixel_size = [1000.00,1000.00]    ENVI_DOIT, 'envi_register_doit', w_fid=fid, w_pos=pos,w_dims=dims,method=7,$     out_name=out_name,pts=ogcp, proj=oproj,r_fid=r_fid3,pixel_size=pixel_size,OUT_BNAME=OUT_BNAME    ENVI_FILE_MNG, id=r_fid3, /remove    ENVI_FILE_MNG, id=fid, /remove  END  

相關推薦

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

【Python學習日記】B站小甲魚：永久儲存（pickle模組）和異常處理（exception）

Pickle 永久儲存 pickle模組用來記錄一些列表，字典，集合等資料，將它們儲存在一個資料夾中作為資料包，等到使用的時候再開啟讀取。