IDL遥感应用入门（13）：采用IDL进行几何校正

这次的文章将介绍一个比较实用的操作，采用IDL进行几何校正操作。几何校正是做遥感图像处理的时候经常涉及到的工作，本文章将系统的讲解整个几何校正的思想和实现步骤。

一些交代

本文的测试数据，链接：http://pan.baidu.com/s/1eQgp0K6 密码：i1ju。关于几何校正，本文只是提到了一个简单的思路，一些特殊的情况，将在最后做一个说明。

读取原始数据

我们可以发现这是一个变形的图像，如果你仔细一点看的话，或许可以发现山东半岛和朝鲜半岛。但是这显示是不符合我们要求的图像，下面说下校正的思路。

校正思路

1. 首先我们需要知道图像中，每一个像素对应的经纬度，这个可以从原来的图像中获取（原来的图像是一景HDF4数据）
2. 求出最大最小的经度值和纬度值，并且求出最大最小经纬度的差值
3. 开辟一个与经纬度差值等比例的数组，然后每一个经纬度做一个映射

校正过程（全部的代码）

;Author:Sailor
;2014-01-15

PRO Course_12

;定义文件路径
MyRootDir='D:\9\'

;遍历文件夹
filearr = file_search(MyRootDir+'TMI','*.hdf',count=num);
FOR fileindex=0,num-1,1 DO BEGIN

;打开SDS模式的HDF文件&返回文件信息
  hdfid=hdf_sd_start(filearr[fileindex])
  hdf_sd_fileinfo,hdfid,datasets,attributes

  print,datasets,attributes
  ;遍历每个属性
  FOR  attrindex=0,attributes-1 DO BEGIN                               
    hdf_sd_attrinfo,hdfid,attrindex,name=a,data=b,count=c,hdf_type=d,type=e
    print,'Name',':',a
    print,'Data',':',b
    print,'Count',':',c
    print,'Hdf_type',':',d
    print,'Type',':',e
    print,'+++++++++++++++++++++++++++++++++++++'
  ENDFOR

  ;遍历每个数据集
  FOR dsindex=0,datasets-1 DO BEGIN

    ;读取数据集数据
    varid=hdf_sd_select(hdfid,dsindex)
    hdf_sd_getdata,varid,Data
    hdf_sd_getinfo,varid,name=dstName,natts=dstNum

    ;如果数据集不为0
    ;IF(dstNum NE 0) THEN BEGIN
      print,'Name',':',dstName
      help,Data
      print,'Num',':',dstNum
      FOR ind2=0,dstNum-1 DO BEGIN
        hdf_sd_attrinfo,varid,ind2,name=disName,data=disData
        print,disName
        print,disData                        
      ENDFOR
      print,'+++++++++++++++++++++++++++++++++++++'
      hdf_sd_endaccess,varid
   ; ENDIF
  ENDFOR

  ;读取highResCh数据集
  ind=hdf_sd_nametoindex(hdfid,'highResCh')
  var_id=hdf_sd_select(hdfid,ind)
  hdf_sd_getdata,var_id,Hrc
  Hrc=Hrc/100+100
  hrc1=REFORM(Hrc[0,*,*])
  hrc2=REFORM(Hrc[1,*,*])

  ;显示原始数据
  tem=image(hrc1,rgb_table=39,title=filearr(fileindex),grid_units=1,POSITION=[0.1,0.1,0.9,1])
  c1 = COLORBAR(TARGET=tem, POSITION=[0.1,0.1,0.9,0.15],ORIENTATION=0,TITLE='bright temperature (K) ')

  ;几何校正
  geodex=hdf_sd_nametoindex(hdfid,'geolocation');得到geo数据集的索引号
  geoid=hdf_sd_select(hdfid,geodex)
  hdf_sd_getdata,geoid,geocon ;获取经纬度信息

  ;求出左下和右上经纬度
  lat=REFORM(geocon[0,*,*])
  lng=REFORM(geocon[1,*,*])
  latmin=min(lat)
  latmax=max(lat)
  LNGMIN=MIN(LNG)
  lngmax=max(lng)
  print,latmin,latmax,lngmin,lngmax

  ;开辟新的数组
  geonew1=FLTARR(320,100)
  geonew2=FLTARR(320,100)
  geonew1[*,*]=0
  geonew2[*,*]=0

  FOR i=0,182,1 DO BEGIN
    FOR j=0,207,1 DO BEGIN
      geonew1[FIX((lng[j,i]-112)/0.1),FIX((lat[j,i]-30)/0.1)]=hrc1[j,i]
      geonew2[FIX((lng[j,i]-112)/0.1),FIX((lat[j,i]-30)/0.1)]=hrc2[j,i]
    ENDFOR
  ENDFOR

  ;过滤数据
  zeidx=where(geonew1 eq 0)
  geonew1[zeidx]=max(geonew1)
  zeidx=where(geonew2 eq 0)
  geonew2[zeidx]=max(geonew2)

  Img=image(geonew1,rgb_table=39,title=filearr(fileindex),grid_units=1,POSITION=[0.1,0.1,0.9,1])
  xaxis=axis('X',location=[0,0],range=[112,143],minor=0, major=17,coord_transform=[112,0.096875],title='Longitude(°)')
  yaxis=axis('Y',location=[0,0],range=[30,39],minor=0, major=7,coord_transform=[30,0.09],title='Latitude(°)')
  c1 = COLORBAR(TARGET=Img, POSITION=[0.1,0.1,0.9,0.15],ORIENTATION=0,TITLE='bright temperature (K) ')

  Img2=image(geonew2,rgb_table=39,title=filearr(fileindex),grid_units=1,POSITION=[0.1,0.1,0.9,1])
  xaxis=axis('X',location=[0,0],range=[112,143],minor=0, major=17,coord_transform=[112,0.096875],title='Longitude(°)')
  yaxis=axis('Y',location=[0,0],range=[30,39],minor=0, major=7,coord_transform=[30,0.09],title='Latitude(°)')
  c2 = COLORBAR(TARGET=Img2,POSITION=[0.1,0.1,0.9,0.15], ORIENTATION=0,TITLE='bright temperature (K) ')

ENDFOR
END

校正结果

一些说明

其实几何校正会涉及很多复杂情况，比如本图中有陆地有海，而在映射的过程中很可能会出现两个点映射到同一个位置的情况，这时候可以做一个平均值，也可以取第一个值或者最后一个值。如果要取平均的话，要注意如果映射到一个点的既有陆地又有海洋，就不应该取平均，这样的均值反而是错误的，遇到这种情况可以查看其邻域的值的大小，如果陆地占多数，就要去掉海洋的值，反之亦然。