IDL对HDF5图像做GLT校正

最近需要对HY-1C的HDF5图像做几何校正，并简单的显示，但是查找了大量的资料，并没有找到相应的代码可供参考，只找到了一位朋友的博客，是对HDF影像做处理，在此感谢，地址在这儿http://blog.sina.com.cn/s/blog_afada1270101b97q.html

由于之前没有接触过IDL开发，于是我一边参考ENVI的帮助文档，一边自己摸索，终于满足需求，代码如下，希望大家做指正。（里面用到了大量的ENVI函数，原生的IDL开发，时间成本有点高，就没有采用）

PRO ENVI_GLT_DOIT_RECORD,_extra=extra
end

pro ENVI_GEOREFERENCE_FROM_GLT_DOIT_RECORD,extra=extra
end

PRO ENVI_GEOREF_FROM_GLT_DOIT_RECORD,extra=extra
end
;上边的代码，是因为我后期用到了C#\idl混合开发，用来处理错误的（参考董彦卿老师的博客）


pro readH5,Fname,WorkSpace
  
  compile_opt idl2
  ;处理ENVI命令，但不显示envi窗体
  e=ENVI(/headless)
  ;/on参数：显示处理进度条
  envi_batch_status_window, /on 
  ;设置输入投影
  inpro=ENVI_PROJ_CREATE(/geographic, datum="WGS-84")
  ;设置输出投影
  outpro=ENVI_PROJ_CREATE(/UTM, zone=51, datum="WGS-84")
  
  ;文件名
;Fname="C:\Users\ASUS\Desktop\1111\H1C_OPER_OCT_L1B_20191111T025000_20191111T025500_06164_10.h5"  

  file_id=H5F_open(Fname)
  ;打开群组
  groupData_id=H5G_OPEN(file_id,'Geophysical Data')
  groupGeo_id=H5G_OPEN(file_id,'Navigation Data')
  ;打开数据集
  datasetR_id=H5D_OPEN(groupData_id,'L_670')
  datasetG_id=H5D_OPEN(groupData_id,'L_865')
  datasetB_id=H5D_OPEN(groupData_id,'L_490')
  
  datasetLon_id=H5D_OPEN(groupGeo_id,'Longitude')
  datasetLat_id=H5D_OPEN(groupGeo_id,'Latitude')
  ;print,datasetLon_id,'  Lon_id'
  ;读取数据集中的shu'ju 
  sdsDataR=H5D_READ(datasetR_id)
  sdsDataG=H5D_READ(datasetG_id)
  sdsDataB=H5D_READ(datasetB_id)
  
  ;[206:1205]这个范围，是我所需的图像处理范围
  sdsDataRChnl=FLOAT(sdsDataR[206:1205,54:1458])
  sdsDataGChnl=FLOAT(sdsDataG[206:1205,54:1458])
  sdsDataBChnl=FLOAT(sdsDataB[206:1205,54:1458])
  sdsDataR=!null
  sdsDataG=!null
  sdsDataB=!null
  
  sdsDataLon=H5D_READ(datasetLon_id)
  sdsDataLat=H5D_READ(datasetLat_id)
  ;Help,sdsDataR
  ;help,sdsDataG
  ;help,sdsDataB
  
  LatChnl=FLOAT(sdsDataLat[206:1205,54:1458])
  LonChnl=FLOAT(sdsDataLon[206:1205,54:1458])
  sdsDataLon=!null
  sdsDataLat=!null
  
  ;处理Latitude波段
  ;判断文件是否存在，若存在，删除重新生成 
  ;RefLatTempFn=e.GetTemporaryFilename()
  RefLatTempFn=WorkSpace+'\RefLatTempFn.dat'
  resLat=file_test(RefLatTempFn)
  
  if resLat eq 1 then begin
    file_delete,RefLatTempFn
    ;file_Dir=file_Dirname(RefLatTempFn)
    RefLatTempHdr=WorkSpace+'\'+'RefLatTempFn.hdr'
    if file_test(RefLatTempHdr) eq 1 then begin
      file_delete,RefLatTempHdr
    endif
  endif
  
  RefLat=ENVIRaster(LatChnl,uri=RefLatTempFn)
  RefLat.Save
  ;print,'*******&&&&&&&&&*******'
  ;help,RefLat
  RefLat_id=ENVIRasterToFID(RefLat)
  LatChnl=!null

  ;处理Longitude波段
  ;RefLonTempFn=e.GetTemporaryFilename()
  RefLonTempFn=WorkSpace+'\RefLonTempFn.dat'
  resLon=file_test(RefLonTempFn)
  if resLon eq 1 then begin
    file_delete,RefLonTempFn
    Lonfile_dir=file_Dirname(RefLonTempFn)
    RefLonTempHdr=WorkSpace+'\'+'RefLonTempFn.hdr'
    if file_test(RefLonTempHdr) eq 1 then begin
      file_delete,RefLonTempHdr
    endif
  endif
  
  RefLon=ENVIRaster(LonChnl,uri=RefLonTempFn)
  RefLon.Save
  RefLon_id=ENVIRasterToFID(RefLon)
  ;print,RefLon_id
  ;help,RefLon_id
  LonChnl=!null
  

  ;生成数据数组，一维：R；二维：G；三维：B
  GeoDataChnl=MAKE_ARRAY([1000,1405,3],TYPE=4)
  FOR i=0,999 DO begin
    FOR j=0,1404 do begin
      ;R
      GeoDataChnl[i,j,0]=sdsDataRChnl[i,j]
      ;G
      GeoDataChnl[i,j,1]=sdsDataGChnl[i,j]
      ;B
      GeoDataChnl[i,j,2]=sdsDataBChnl[i,j]
    endfor
  endfor
 ;help,GeoDataChnl
 
 ;RefGeoDataTempFn=e.GetTemporaryFilename()
 RefGeoDataTempFn=WorkSpace+'\RefGeodata.dat'
 resGeoData=file_test(RefGeoDataTempFn)
 if resGeoData eq 1 then begin
  file_delete,RefGeoDataTempFn
  GeoDataFile_Dir=file_Dirname(RefGeoDataTempFn)
  RefGeoDataTempHdr=WorkSpace+'\'+'RefGeodata.hdr'
  if file_test(RefGeoDataTempHdr) eq 1 then begin
    file_delete,RefGeoDataTempHdr
  endif
 endif
    
 RefGeoData=ENVIRaster(GeoDataChnl,uri=RefGeoDataTempFn,interleave='BSQ')
 RefGeoData.Save
 RefGeoData_id=ENVIRasterToFID(RefGeoData)
 ;print,RefGeoData
 ;print,RefGeoData_id,'   RefGeoData_id'
 ;获取refGeoData_dims，用于下一步做GLT
 ENVI_FILE_QUERY,RefGeoData_id,dims=refGeoData_dims
 ;print,'*********overRead*************' 
 GeoDataChnl=!null
 
 ;生成GLT
 ;gltTempFn=e.GetTemporaryFileName()
 gltTempFn=WorkSpace+'\glt.dat'  
 ;RUN SOME APLICATION
 ENVI_DOIT,'ENVI_GLT_DOIT',DIMS=refGeoData_dims,I_PROJ=inpro,O_PROJ=outpro,$
    OUTNAME=gltTempFn,R_FID=glt_id,ROTATION=0,X_FID=RefLon_id,Y_FID=RefLat_id,$
    X_POS=[0],Y_POS=[0]
 ;print,'*********overGLT*************'
 
 
 glt_outname=WorkSpace+'\END_GLT.dat'
 resEndGlt=file_test(glt_outname)
 if resEndGlt eq 1 then begin
  file_delete,glt_outname
  ;endGltDir=file_Dirname(glt_outname)
  glt_outHdr=WorkSpace+'\'+'END_GLT.hdr'
  if file_test(glt_outHdr) eq 1 then begin
    file_delete,glt_outHdr
  endif
 endif
 

  ;做glt校正
 ENVI_DOIT,'ENVI_GEOREF_FROM_GLT_DOIT',$
      BACKGROUND=0,FID=RefGeoData_id,GLT_FID=glt_id,OUT_NAME=glt_outname,POS=[0:2]
      
 ;print,glt_outname
 ;print,'*********over*************'

  ;void=DIALOG_MESSAGE(!ERROR_STATE.MSG,/INFOR,TITLE='CW')
end

参考：http://www.harrisgeospatial.com/docs/ENVI_FILE_QUERY.html#DATA_TYPE

因为之前做AE开发的时候，就是有的功能不会做，就去用ArcGis desktop把思路走通，然后再用Arcengine做。这次做GLT几何校正，没想到也能实现。一些小感想：idl开发感觉就和matlab差不多，语言风格都相似，然而发现用idl做研究的居然不多，一位大佬告诉我，说国内外普遍用matlab+C++实现，要不这么难找idl资料呢。

下一篇文章，简单阐述一下C#+idl混合开发