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