ENVI/IDL——擷取AVIRIS資料波段值和經緯度

以前，這一步，我都是手動導出的，詳細請看另一篇博文：http://blog.chinaunix.net/uid-28490468-id-3462577.html

最近要處理較多的遙感資料，是以要把這一步通過代碼自動化。

看了看help文檔，發現ENVI支援一種叫IDL的語言，完全沒有聽說過，可能是我孤陋寡聞了，呵呵，有同感的舉手。

有好玩的東西，那麼就開始學呗。

打開ENVI+IDLX.X，編寫IDL的IDE會自動和ENVI一起打開，這IDE很面熟，像eclipse。

IDL有自己的help文檔，通過IDE可以打開。

另外網上也有些資料：

簡潔明了： http://www2.geog.ucl.ac.uk/~mdisney/teaching/unix/idl/idl.html 

一個國内的論壇，要注冊，稽核時間要幾天（最好快點，呵呵）：http://bbs.esrichina-bj.cn/ESRI/forumdisplay.php?fid=28

各種資料胡亂看了一通，稍微總結幾點：

IDL

1）IDL的數組是列優先的array[cols, rows]

2）數組是0開始的，和matlab不一樣

3）procedure和function都叫routine，其中檔案命名

File naming—a program file must be named the same as the main routine.

4）Primary routine must be the last in the file，這個關系到自動編譯

5）對routine的參數傳遞有兩種，argument和keyword

具體大家看help文檔吧。

但是光有IDL這些基本的知識，根本無法處理像AVIRIS這種超光譜圖檔，我甚至嘗試着去手動解析hdr頭檔案。

不過在愚蠢之餘，試着google大法了一把，結果發現原來ENVI用IDL另外開發了一些API，特意用來處理碩大的遙感圖檔，

太激動了，呵呵。http://www.360doc.com/content/10/1101/20/472115_65790354.shtml

這些API在IDL的help文檔裡沒有，要在ENVI的help文檔裡有，真是坑爹。

ENVI+IDL

怎麼使用這些API來完成我的任務呢？我是參考：http://www.360doc.com/content/10/0405/21/472115_21741788.shtml

用第一題的做法就能完成我的任務了，不過讀取波段值可以用另一個函數，更友善。

具體API

擷取波段值 envi_open_file      filename, r_fid=fid, /no_realize, /no_interactive_query

     if (fid eq - 1) then return

envi_file_query      fid, ns=ns, nl=nl, nb=nb, wl=wl

envi_get_slice

擷取經緯度 envi_convert_file_coordinates envi_convert_projection_coordinates envi_get_projection envi_proj_create

下面是我的實驗代碼：（要把圖檔和對應的hdr頭檔案放在一個檔案目錄下）

PRO file_test
  data_dir = "E:\Data\OilLeakData\data\f100524t01p00r11rdn_b\"
  filename = data_dir + "f100524t01p00r11rdn_b_sc01_ort_img"
  envi_open_file, filename, r_fid=fid, /no_realize, /no_interactive_query
  if (fid eq -1) then return ;;check fid
  envi_file_query, fid, ns=ns, nl=nl, nb=nb
  pos = lindgen(nb);
  xs = 225
  xe = 344
  ys = 6183
  ye = 6193
  iproj = envi_get_projection(fid=fid)
  oproj = envi_proj_create(/geographic)
  print, 'content'
  for y = ys, ye do begin
    for x = xs, xe do begin
      data = envi_get_slice(fid=fid, line=y, pos=pos, $
        xs=x, xe=x, /bip)   ;; get bands' value
      envi_convert_file_coordinates, fid, x, y, xmap, ymap, /to_map
      envi_convert_projection_coordinates, xmap, ymap, iproj, oxmap, oymap, oproj ;;Lon = oxmap, Lat = oymap
      help, x, xmap, oxmap, data
      print, x, y, xmap, ymap, oxmap, oymap, data
    endfor
  endfor
  envi_file_mng, id=fid, /remove ;;Don't forget to close file
END