IDL考前复习(五)IDL联动ENVI辐射定标

前言

  千呼万唤始出来，一直运行不成功原来是图像的问题，打不开第二波段的图像，人都给你整傻咯，IDL是真心的不好用而且还要和ENVI联动，就更蠢了，ENVI自有的辐射定标功能它不香吗，可是老师上课写了考试还考，只能含泪写完了它,流下了技术人debug一天的眼泪。

需求分析

虽然流程看起来比较简单，但是中间的文本的操作真的是复杂到哭，真的是含泪写完了博客。

1.打开txt文件利用hash表存储所要的信息
  何是 $hash，hash$ 表就是一个键对于一个值，在下面的部分我们用 $hsmeta[linekey]=linevalue$ 来表示。由于读取中有很多的复杂的操作，我尽量也讲清楚一些，用读取出的一行显示的截图来说明一些问题。

function  getmetadata,filename
openr,lun,filename,/get_lun             ;
line=''
hsmeta=hash()
  while(~eof(lun)) do begin
  readf,lun,line
  if(line.indexof('=')ge 0) then begin
      linearray=line.SPLIT('=')
      linekey=linearray[0].trim()
      linevalue=linearray[1].trim()
      hsmeta[linekey]=linevalue
   endif
  endwhile
return,hsmeta
end

要得到所需要的键和值的对应关系，打开文件读文本的操作我也就不解释了，看过以前的blog应该就明白，然后 $line$ 是我定义好的文本变量，接收我逐行读取的文本， $hsmeta$ 是一个hash表，后面开始循环逐行读文本， $readf,lun,line$ 应该都很熟悉了。这个函数输入参数为文件名，输出的是hsmeta一个大的 $hash$表

   if(line.indexof(’=’)ge 0) then begin

为什么要设置这个循环呢，是因为那个MTL.TXT文件里最后一行是以 end 结尾的，如果读到那里会导致一些错误，即我希望我的hash表里都是下面这样的对应关系

$$REFLECTANCE\_ADD\_BAND\_1=-0.100000$$

而如果不加入，hash表中最后会出现一个end，会出现键和值不对应的问题。

   end
循环是个整体框架，循环中带有一些文本的操作，即：

  linearray=line.SPLIT(’=’)
   linekey=linearray[0].trim()
   linevalue=linearray[1].trim()

我怕解释不清楚就以截图来说明，首先用split函数把line切割开来，把 $linearray$ 切割成2个字符串部分,看图片很清晰，而后面 $linekey=linearray[0].trim$ 首先我取出了第一个元素，而这个元素会有很多空格，我需要用 $trim$ 函数把空格全部删掉，同理 $linevalue$ 也是一样的。

下面来最麻烦的东西了就是 $hash$,$hash$ 一个键对应一个值，例如下面

$hash[linekey]$ 是一个键，而对应的是 $linevalue$ 这么一个值，我们可以通过一个键对应一个值，所以后续其实就是无数次循环，把整个存成一个大的 $hash$ 以便后续的操作。最后这个函数返回了 $hsmeta$ 以便后续的函数调用。

2.联动ENVI进行相关的操作

  其实本来是很简单的一件事情，只不过正好在学IDL和ENVI混合编程，所以老师一定要用ENVI进行操作也没办法，那就肝呗，其实核心就是
$$ra{d}_{DN}=DN\ast gain+bais$$罢了，有点废话了干活。

function  getcalculated,band,hashmeta,folder
  bandstr=string(band)
  bandstr=bandstr.trim()
  gainkey='REFLECTANCE_MULT_BAND_'+bandstr
  biaskey='REFLECTANCE_ADD_BAND_'+bandstr
  filename='FILE_NAME_BAND_'+bandstr
  
  gains=float(hashmeta[gainkey])
  bias=float(hashmeta[biaskey])
  fn=hashmeta[filename]
  fn=folder+fn.replace('"','')
  
  envi_open_data_file,fn,r_fid=fid,/tiff
  envi_file_query,fid,dims=dims,map_info=mapinfo,nl=nl,ns=ns,nb=nb,bnames=bnames,interleave=interleave,offset=offset
  
  mapinfo=envi_get_map_info(fid=fid)
 
  lsdata=envi_get_data(fid=fid,dims=dims,pos=0)
  result=lsdata*gains+bias
  envi_write_envi_file,result,out_name='E:/radcalib.img',map_info=mapinfo,data_type=4,nl=nl,ns=ns,nb=nb,bnames=bnames,interleave=interleave,offset=offset
return,'ok'
end

这个函数的传入值是 $band$ ，即你要处理的波段。

   bandstr=string(band)
   bandstr=bandstr.trim()

看图说话还是比较容易的, $band$ 在实际函数传递的时候是个int类型的,所以我需要把它转成 string 的文本类型而转化成文本，而前面又有空格所以就需要用 trim 进行裁剪，

   gainkey=‘REFLECTANCE_MULT_BAND_’+bandstr
   biaskey=‘REFLECTANCE_ADD_BAND_’+bandstr
   filename=‘FILE_NAME_BAND_’+bandstr

这段我就不大想解释了，其实就是一些字符串的命名 ‘+’ 号就是把字符串给拼接起来，我就不截图了，太伤身体了。

   gains=float(hashmeta[gainkey])
   bias=float(hashmeta[biaskey])

由于在之前那个函数中已经有一个 $hashmeta$ 了，其实现在都是一些形式参数，等到了主pro过程是要传入实际参数的。
在主pro过程中会先利用第一个函数 $hashmeta=getmetadata(folder+filename)$ 会得到一些键和一些值。
就是说找到相应的键和值由于那些键和值都是文本类型的，需要转化成浮点型的。

   fn=hashmeta[filename]
   fn=folder+fn.replace(’ " ‘,’’)

这里就是找到文件名，根据 $hash$表 找到文件名，然后把双引号替代成单引号，并且加上它的文件夹，我特意断点截图了一下，看了一下效果。

中途的过程太过心酸，实际上单纯看这个函数定义没啥用，因为没传入实参，等到主pro过程加入实参会更有用，其实就是根据 $hash$表 相应的键找到相应的值。看到这个过程的原因是我已经开始在主pro过程运行完第一个函数了。

envi_open_data_file,fn,r_fid=fid,/tiff
envi_file_query,fid,dims=dims,map_info=mapinfo,nl=nl,ns=ns,nb=nb,bnames=bnames,interleave=interleave,offset=offset

这部分markdown没法一行用茶青色和特定的字体做完只好贴代码了。

开始了envi联动的过程，记得如果envi联动的话需要在命令行输入ENVI，打开ENVI经典的进程，否则一定运行不出来

解释下，就是打开利用envi打开文件，文件名，然后记得是 r_fid = fid 后面有好多 fid 或者 fid = fid 的，并且打开的是tiff文件。

下一个就是利用ENVI查询相应的参数，这么多参数我实在没法解释了，我猜也不能考把，知道mapinfo，interleave 和 offset 应该就够了，投影、乘常数、偏移量。

   mapinfo=envi_get_map_info(fid=fid)
   lsdata=envi_get_data(fid=fid,dims=dims,pos=0)

第一句话就是得到投影的意思，上面的 query 只能查询好像不能得到参数，需要专门的 envi_get_map_info 函数得到图像的地理参考信息，而且记得是 fid = fid 上面用查询 query 函数的时候只要一个 fid 不需要额外的添加。

第二句话就是得到数据的意思， fid = fid 注意下， 因为只有一个波段所以设置 pos = 0。

   result=lsdata*gains+bias

这句话才是整个程序的核心，就是辐射定标的公式，也就是一个矩阵计算罢了。

envi_write_envi_file,result,out_name='E:/radcalib.img',map_info=mapinfo,data_type=4,nl=nl,ns=ns,nb=nb,bnames=bnames,interleave=interleave,offset=offset

最后写出ENVI文件，记得写出的数据的值是 float 类型即 datatype = 4，然后又是一堆烦躁的参数。

   return,‘ok’
这个我觉得是老师做了10多年的程序主管习惯了，就只是告诉你程序运行成功了。

3.主pro过程

终于来到最后了，我觉得不仅仅是我最后用IDL了，也可能是我最后一段时间再做普通遥感的内容了，毕设以及后期应该会转图像处理和机器学习的内容，伤感一下，哈哈，算半逃出升天把进入另一个魔窟吧。

pro  fushedingbiao
folder= 'E:\LC81220352018123LGN00\'
filename='LC08_L1TP_122035_20180417_20180501_01_T1_MTL.txt'
hashmeta=getmetadata(folder+filename)
result=getcalculated(3,hashmeta,folder)
print,result
end

过程的名称，文件夹的名称、文件夹、调用两个函数，打印结果。

看看吧，运行成功就会有这些东西。

写在最后的话

  IDL系列应该是更新完了，帮助组内的胖友过考试应该没啥问题了，没啥要求，有杯肥宅快乐水就好了。
  开玩笑，大四了，距离考研时间也不多了，但是作为一名技术狂热粉，还是一个强迫症患者，不运行成功不写出blog实在太难受了，所以就花了下午2个小时写完了这篇blog，对了想起来了,IDL的三维彩色图像显示是个ZZ，列和波段(MATLAB说页)会偷偷地互换，至于咋互换，慢慢探索去吧，哈哈哈~~！！
  写博客真的可以放飞自我，如果没人发现我的话，可以把自己很多想说的话都留在博客中，我不喜欢当一个 institutionalized 的人，大学本身就应该是一个自我探索兴趣的过程，如果还把自己禁锢在背书，不实践探索自己兴趣的时段，那么应该人生不会再有机会了把。取得好成绩和找到好兴趣可能可以兼得，反正我应该没能兼得，但是如果你问我后悔吗，不后悔，人生要是能后悔的话，那么人人都应该在清华北大了吧，最后的最后，IDL再见，资源遥感、应用遥感再见！