基于IDL的高分二号影像批量预处理程序(第二版)

上一版本插件预处理效果不太理想(特别是融合后分辨率不符合预期),上一版资源不再提供下载,但处理思路可以借鉴。目前在本机40GB RAM工作站上预处理一景影像大概需要13分钟。

废话少说先上资源
如果各位看官支持原创、积分充盈,可以点击下方链接使用积分下载资源
https://download.csdn.net/download/weixin_43490758/85794749

想白嫖也没关系,百度云资源也双手奉上,链接及提取码见页面底端

提示: 如果在运行过程中出现个别问题,建议先认真浏览正文,一般情况下能够解决绝大多数问题,实在解决不了的再在评论区提问,当个社畜着实不易,不能做到及时回复还请见谅!

说点真心话
目前还有很多功能需要完善,我也是IDL新手,需要参考网上相关资源,功能也在不断完善和增加,后期工作不太忙的话可能会研究批量空间配准和镶嵌匀色等功能,敬请期待。目前手头积分余额不足,所以厚着脸皮来赚点积分方便后期功能拓展。保证下面所提及的所有功能均已实现,亲测有效。

一、程序界面及主要功能

目前本程序仅针对国产高分二号遥感影像的批预处理,导出数据格式为dat。按照处理流程依次为①文件解压(可选)②辐射定标(全色、多光谱)③快速大气校正(多光谱)④几何校正(全色、多光谱)⑤融合。 程序功能界面如下图所示:
插件主界面

输入文件夹】应包含所有需要预处理的原始高分二号影像(扩展名为*.tar.gz或*.tiff),建议不要对下载的影像做任何编辑。如果影像未解压,程序会在输入文件夹下创建同名文件夹并将文件解压至该文件夹内。
国产插件sav】ENVI APP store下载的国产卫星支持插件,全名为(envi_china_satellites_support.sav),在应用商店自行下载或在本人上传的资源中下载,选择sav文件即可。貌似该工具只支持5.3以上版本,低版本用户请谨慎使用。
保留近红外波段】可以通过单选按钮选择是否保留近红外波段(舍弃近红外波段可以在一定程度上提高处理效率,之前做过测试一景影像大概能缩减10分钟时间),默认保留近红外波段。
融合方式】目前只提供了NNDiffusePanSharpening、GramSchmidtPanSharpening两种融合方式,默认为NNDiffusePanSharpening。
输出文件夹】所有过程文件(文件名以temp结尾)及最终结果文件会存储在该文件夹中,当一景影像处理完成后程序会自动将过程文件删除。(大批量的遥感影像占用内存大,建议将该文件夹放在具有足够内存空间的磁盘)

二、处理结果展示

话不多说,直接上图看效果。
图:原始全色影像(1m分辨率)
原始全色影像
图:原始多光谱影像(4m分辨率)
原始多光谱影像
图:预处理后融合影像(1m分辨率)
批量处理后融合影像

三、程序环境要求

目前使用的软件是ENVI5.3(64-bit),IDL8.5,win10操作系统,win11操作系统中会出现程序界面控件显示不全。并不能保证在其他环境下能够稳定有效运行。

四、功能模块代码

①程序主入口、界面

;高分影像批量预处理
PRO PreTreat
  defsysv,'!mosicmethod','NNDiffusePanSharpening'
  defsysv,'!RGBN',1
  COMPILE_OPT IDL2
  
  ;为IDL编译器指定编译规则,idl2是(DEFINT32,STRICTARR)的简写,DEFINT32为IDL指定整型为32位,而不是16位。STRICTARR则强制IDL编译器以[]作为数组的定义。
  ;
  ;
  ;/headless  不显示envi界面

  ;创建用户界面
  ;
  ;主窗体
  wtop = widget_base(title='国产高分影像批量预处理(V2.0)', xsize=300, ysize=300, /column, /TLB_KILL_REQUEST_EVENTS,TLB_FRAME_ATTR=1, mbar=menubar)
  ; /column 表示列数为1,  多行每行左右顶到头      /row  多列
  ; TLB_KILL_REQUEST_EVENTS  是否返回关闭事件
  ; TLB_FRAME_ATTR=1创建的窗口类型,1表示窗口无法进行大小、最大化操作
  ; mbar=menubar  传出菜单栏的ID
  
  ;选择影像所在文件夹
  winputBase = WIDGET_BASE(wtop,  XPAD=10, YPAD=5,  /FRAME, /row)
  ;/FRAME 部件边框宽度
  winputLabel = WIDGET_LABEL(winputBase, VALUE='输入文件夹:')
  winputTextField = widget_text(winputBase, EDITABLE=0 ,UVALUE='InputTEXT')
  winputbutton=widget_button(winputBase, VALUE='选择',UVALUE='InputButton')
  
  ;选择国产高分影像支持插件sav文件
  wChinasupportBase = WIDGET_BASE(wtop,  XPAD=10, YPAD=5,  /FRAME, /row)
  wChinasupportLabel = WIDGET_LABEL(wChinasupportBase, VALUE='国产插件sav:')
  wChinasupportTextField = widget_text(wChinasupportBase, EDITABLE=0 ,UVALUE='ChinasupportTEXT')
  wChinasupportbutton=widget_button(wChinasupportBase, VALUE='选择',UVALUE='ChinasupportButton')

  ;选择处理波段数
  wcheckframe = WIDGET_BASE(wtop, ysize=40,/BASE_ALIGN_CENTER, XPAD=10, YPAD=5, /FRAME, /row)
  wcheckLabel = WIDGET_LABEL(wcheckframe, /align_center, VALUE='保留近红外波段:')
  wcheckBase = WIDGET_BASE(wcheckframe,/BASE_ALIGN_CENTER, XPAD=10, YPAD=5,/EXCLUSIVE,/ROW)
  RGB = WIDGET_BUTTON(wcheckBase,value ='是',uName = 'checkright',UVALUE='RGBN', /NO_RELEASE )
  RGBN = WIDGET_BUTTON(wcheckBase,value ='否',uName = 'checkerror',UVALUE='RGB',/NO_RELEASE )
  WIDGET_CONTROL, RGB, /SET_BUTTOn

  ;选择融合方式
  wmosicbase = WIDGET_BASE(wtop, /BASE_ALIGN_CENTER, XPAD=10, YPAD=5, /FRAME, /row)
  wmosicLabel = WIDGET_LABEL(wmosicbase, /align_center,VALUE='选择融合方式:')
  wmosicDroplist = WIDGET_DROPLIST(wmosicbase, VALUE=['NNDiffusePanSharpening', 'GramSchmidtPanSharpening'], UVALUE='MOSICLIST')
  WIDGET_CONTROL, wmosicDroplist, SET_DROPLIST_SELECT=0;默认融合方式为NNDiffusePanSharpening
  
  ;选择影像输出文件夹
  woutputBase = WIDGET_BASE(wtop, /BASE_ALIGN_CENTER, XPAD=10, YPAD=5, /FRAME, /row)
  woutputLabel = WIDGET_LABEL(woutputBase, VALUE='输出文件夹:')
  woutputTextField = widget_text(woutputBase, EDITABLE=0, UVALUE='OutputTEXT')
  woutputbutton=widget_button(woutputBase, VALUE='选择', UVALUE='OutputButton')

  ;执行命令
  wbuttonBase = WIDGET_BASE(wtop, /BASE_ALIGN_BOTTOM, XPAD=90, YPAD=10, /FRAME, /row)
  wsurebutton=widget_button(wbuttonBase, /ALIGN_CENTER, xsize=50, ysize=25,UVALUE='okButton',VALUE='确定')
  wcancelbutton=widget_button(wbuttonBase, /ALIGN_CENTER,xsize=50, ysize=25,UVALUE='cancelButton',VALUE='关闭')
  
  wtipBase = WIDGET_BASE(wtop, /BASE_ALIGN_BOTTOM, XPAD=10, YPAD=5, /FRAME, /row)
  wtipLabel = WIDGET_LABEL(wtipBase, VALUE='提示:请务必使用ENVI5.3以上版本软件!!!')
  ; 实例化控制顶级窗体及其子控件的显示
  WIDGET_CONTROL, wtop, /REALIZE
  
 ;创建结构体,包含各个组件ID 和参数
  pState={INPUTTXT:winputTextField,$
    ChinasupportTXT:wChinasupportTextField,$
    OUTPUTTXT:woutputTextField,$
    MOSICLIST:wmosicDroplist}
    
  WIDGET_CONTROL, wtop,set_uvalue=ptr_new(pState)
  ;绑定事件
  Xmanager,'PreTreat',wtop, /no_block
END

②响应事件

;事件响应
pro PreTreat_Event, Event     ; IN: event structure
  e = ENVI(/headless)
  ;获取top指针
  
  WIDGET_CONTROL, Event.top, GET_UVALUE=sState
  ;  界面右上角关闭按钮可用
  if (TAG_NAMES(Event, /STRUCTURE_NAME) EQ $
    'WIDGET_KILL_REQUEST') then begin
    result=DIALOG_MESSAGE("是否关闭",/QUESTION, TITLE='提示')
    CASE (result) OF
      'Yes': BEGIN
        PTR_FREE, sState;销毁指针
        WIDGET_CONTROL,event.TOP,/DESTROY
        return
      END
      'No': BEGIN
        RETURN
      END
      ELSE: BEGIN
      END
    ENDCASE
  endif

 
  ;根据点击对象的ID值获取 其UVALUE值
  WIDGET_CONTROL, Event.id, GET_UVALUE=eventval
  case eventval of
    'InputButton' : begin 
      inFile = DIALOG_PICKFILE(DIALOG_PARENT = e.WIDGET_ID, /DIRECTORY, TITLE='选择原始影像所在文件夹')
      ;如果点击了取消或者文件无效,跳出
      if ~file_test(inFile) then return
      ;widget_control, sState.wTextField,  get_value=textVal
      WIDGET_CONTROL,(*sState).INPUTTXT, set_Value = inFile
      print,inFile
    end ;   

    'ChinasupportButton' : begin
      ChinasupportFile = DIALOG_PICKFILE(DIALOG_PARENT = e.WIDGET_ID, TITLE='选择国产卫星支持sav插件',FILTER='*.sav')
      ;如果点击了取消或者文件无效,跳出
      if ~file_test(ChinasupportFile) then return
      WIDGET_CONTROL,(*sState).ChinasupportTXT, set_Value = ChinasupportFile
      print,ChinasupportFile
    end ;

    'RGB' : begin 
      !RGBN=0
      print,!RGBN
    end

    'RGBN' : begin   
      !RGBN=1
      print,!RGBN
    end
    
    'MOSICLIST' : begin
      ;获取融合方法列表中的所有值存入mosicmethods
      widget_control, (*sState).MOSICLIST,  get_value=mosicmethods
      ;获取下列表选择的索引值
      listValue = WIDGET_INFO(Event.id, /DROPLIST_SELECT)

      ;获取选择的融合方法名
      !mosicmethod=mosicmethods[listValue]
      print,!mosicmethod
    end     
    
    
    'OutputButton' : begin
      outFile = DIALOG_PICKFILE(DIALOG_PARENT = e.WIDGET_ID, TITLE = '选择输出目录', /DIRECTORY)
      if ~file_test(outFile) then return
      WIDGET_CONTROL,(*sState).OUTPUTTXT, set_Value = outFile
      print,outFile
    end ;
    
    'okButton': begin
      WIDGET_CONTROL, (*sState).INPUTTXT, get_value=inputdir
      WIDGET_CONTROL, (*sState).ChinasupportTXT, get_value=Chinasupportfile
      WIDGET_CONTROL, (*sState).OUTPUTTXT, get_value = outputdir
      ;如果输入输出文件均设置完成,就开始执行
      IF inputdir NE "" && outputdir NE "" THEN BEGIN 
        ;检索文件夹中扩展名为tiff的文件,获取数量    num和满足要求的文件名列表filelists
        filelists=FILE_SEARCH(inputdir,'*.tiff',count = num)
        ;这里之所以检索后缀为tiff的文件而不是直接检索压缩包文件是由于,用户可能已经自行完成解压但是没有将压缩包文件删除,否则可能导致多次解压产生重复数据
        ;如果不存在后缀为tiff的文件,先进行解压操作
        IF num eq 0 THEN BEGIN
          filelists=file_search(inputDir,'*.gz',count = num)
          FOR i=0, num-1 DO BEGIN
            print, filelists[i]
            tarout_folder = Untar(filelists[i])
            print, tarout_folder
          ENDFOR
          endfile = PrecessingBatch(inputdir, Chinasupportfile, !mosicmethod, !RGBN, outputdir)
          ;如果有满足条件的文件
        ENDIF ELSE BEGIN
          endfile = PrecessingBatch(inputdir, Chinasupportfile, !mosicmethod, !RGBN, outputdir)
        ENDELSE
        
        
      ;否则弹出警告框进行提示
      ENDIF ELSE BEGIN
        result=DIALOG_MESSAGE("输入或输出路径不能为空", TITLE='参数异常')
        CASE (result) OF
          'OK': BEGIN
            return
          END
          ELSE: BEGIN
          END
        ENDCASE
      ENDELSE
    END
    
    'cancelButton': begin
      result=DIALOG_MESSAGE("是否关闭",/QUESTION, TITLE='提示')
      CASE (result) OF
        'Yes': BEGIN
          PTR_FREE, sState;销毁指针
          WIDGET_CONTROL,event.TOP,/DESTROY
          return
        END
        'No': BEGIN
          return
        END
        ELSE: BEGIN
        END
      ENDCASE
     END
    
     ELSE: begin
      PRINT, '未找到匹配项'
     END
  endcase
end

③单个tar.gz文件解压

;参数说明
;inputFile:单个压缩文件路径;
;tarout_folder:解压目标文件夹;

Function Untar, inputFile, tarout_folder
  tarout_folder = FILE_DIRNAME(inputFile)+'\'+FILE_BASENAME(inputFile,'.tar.gz')
  ;解压GZ,解压后的后缀文件是.tar
  FILE_GUNZIP, inputFile
  ;创建解压后gz同名文件夹,否则默认解压到当前文件夹
  IF ~FILE_TEST(tarout_folder) THEN FILE_MKDIR, tarout_folder
  ;解压tar
  FILE_UNTAR, FILE_DIRNAME(inputFile)+'\'+FILE_BASENAME(inputfile,'.tar.gz')+'.tar',tarout_folder
  ;删除中间解压tar文件
  FILE_DELETE,FILE_DIRNAME(inputFile)+'\'+FILE_BASENAME(inputfile,'.tar.gz')+'.tar'
  ;返回解压后文件夹名
  return,tarout_folder
end

④预处理单景遥感影像

;参数说明
;msfile:多光谱影像路径;
;pnfile:全色影像路径;
;Chinasupportfile:国产卫星支持插件sav文件
;mosicmethod:融合方法(可选值有'NNDiffusePanSharpening''GramSchmidtPanSharpening');
;isrgbn:是否包含近红外波段(可选值有10);
;outputdir:输出文件夹;
;tarout_file:结果文件路径。

Function PrecessingSingle, msfile, pnfile, Chinasupportfile, mosicmethod, isrgbn, outputdir,tarout_file
  COMPILE_OPT IDL2
  e = ENVI(/headless)
  Restore, Chinasupportfile ;需要sav文件完整路径
  
  t = SYSTIME(1)
  tarout_name=strmid(FILE_BASENAME(msfile),31,13)

  ;1_RadiometricCalibration
  ms = ENVI_Open_GF2_Raster(msfile);使用APP store下载的支持国产卫星元数据读取的软件打开  
  rad = ENVITask('RadiometricCalibration') ;ENVI5.2.1 Introduced
  rad.Input_Raster = ms
  rad.CALIBRATION_TYPE='Radiance';
  rad.SCALE_FACTOR=0.1
  rad.OUTPUT_DATA_TYPE='Float'
  rad.Output_Raster_URI = outputdir+'GainOffset_temp.dat'
  rad.EXECUTE

  
  ;pn = e.Openraster(pnfile)
  pn = ENVI_Open_GF2_Raster(pnfile);使用APP store下载的支持国产卫星元数据读取的软件打开  
  
  radpn = ENVITask('RadiometricCalibration')
  radpn.Input_Raster = pn
  radpn.CALIBRATION_TYPE = 'Top-of-Atmosphere Reflectance'   ;'Radiance'; 默认     ;
  radpn.Output_Data_Type = 'UInt'
  radpn.SCALE_FACTOR = 10000
  radpn.Output_Raster_URI = outputdir+'panGainOffset_temp.dat'
  radpn.EXECUTE
  t1=SYSTIME(1)
  PRINT,'完成第一步【辐射定标】 ,耗时' + STRING((t1-t)/60, format='(f7.2)') + '分钟'
  PRINT,'*********************************************************************'
  ;2_QUAC
  PRINT,'正在执行第二步【快速大气校正】'
  ms = rad.Output_Raster
  qac = ENVITask('QUAC') ;ENVI5.1 Introduced
  qac.Input_Raster = ms
  qac.Output_Raster_URI = outputdir+'QUAC_temp.dat'
  qac.EXECUTE
  t2=SYSTIME(1)
  PRINT,'完成第二步【快速大气校正】 ,耗时'$
    + STRING((t2-t1)/60, format='(f7.2)') + '分钟,共耗时' + STRING((t2-t)/60, format='(f7.2)') + '分钟'
  PRINT,'*********************************************************************'

  ;3_RPCOrthorectification
  PRINT,'正在执行第三步【正射校正】'
  
  dmfile = e.Root_Dir + 'data\GMTED2010.jp2'
  dm = e.Openraster(dmfile)
  ort = ENVITask('RPCOrthorectification') ;ENVI5.1 Introduced
  ort.dem_raster = dm
  ms = qac.Output_Raster
  IF isrgbn THEN BEGIN
    ort.input_raster = ms
  ENDIF ELSE BEGIN 
    RGBRaster=ENVISubsetRaster(ms,Bands=[0,1,2]) ;
    ort.input_raster = RGBRaster
  ENDELSE
  ort.Output_Pixel_Size = 4
  ort.Output_Raster_URI = outputdir+'RPCOrthorectification_temp.dat'
  ort.EXECUTE

 
  pn=radpn.Output_Raster
  ortpn = ENVITask('RPCOrthorectification')
  ortpn.dem_raster = dm
  ortpn.input_raster = pn
  ortpn.Output_Pixel_Size = 1
  ortpn.Output_Raster_URI = outputdir+'panRPCOrthorectification_temp.dat'
  ortpn.EXECUTE
  t3=SYSTIME(1)
  PRINT,'完成第三步【正射校正】 ,耗时'$
    + STRING((t3-t2)/60, format='(f7.2)') + '分钟,共耗时' + STRING((t3-t)/60, format='(f7.2)') + '分钟'
  dm.close
  PRINT,'*********************************************************************'


  ;5_GramSchmidtPanSharpening
  PRINT,'正在执行第四步【融合】'
  print, mosicmethod
  ms = ort.Output_Raster
  pn = ortpn.Output_Raster
  gsf = ENVITask(mosicmethod) ;ENVI5.2 Introduced
  gsf.Input_Low_Resolution_Raster = ms
  gsf.Input_High_Resolution_Raster = pn
  tarout_file= outputdir+ tarout_name+'.dat'
  gsf.Output_Raster_URI = tarout_file
  ;gsf.IGNORE_VALIDATE = 1 
  ;gsf.PIXEL_SIZE_RATIO = 4
  gsf.EXECUTE
  t5=systime(1)
  PRINT,'完成第四步【融合】 ,耗时'$
    + STRING((t5-t3)/60, format='(f7.2)') + '分钟,共耗时' + STRING((t5-t)/60, format='(f7.2)') + '分钟'
  PRINT,'*********************************************************************'

  ms.close, error=err
  pn.close, error=err
  rad.OUTPUT_RASTER.Close, error=err
  radpn.OUTPUT_RASTER.Close, error=err
  qac.OUTPUT_RASTER.Close, error=err
  ;rgs.Output_Raster.Close,error=err
  gsf.Output_Raster.Close, error=err
  ort.OUTPUT_RASTER.Close, error=err
  ortpn.OUTPUT_RASTER.Close, error=err
  
  ;循环释放内存
  opendata=e.getopendata()
  length=size(opendata,/dimensions)
  for i=0,length[0]-1 do begin
    opendata[i].close
  endfor
  
 

  ;删除中间临时文件
  filelists=FILE_SEARCH(outputdir,'*temp*',count = num)
  FOR i=0, num-1 DO BEGIN
    FILE_DELETE, string(filelists[i])
  ENDFOR
  RETURN,tarout_file
  e.Close
END

⑤批量预处理影像

;参数说明
;inputdir:待处理影像所在文件夹;
;Chinasupportfile:国产卫星支持插件sav文件;
;mosicmethod:融合方法(可选值有'NNDiffusePanSharpening''GramSchmidtPanSharpening');
;isrgbn:是否包含近红外波段(可选值有10);
;outputdir:输出文件夹
;tarout_file:输出文件路径(只返回最后一景影像路径)

Function PrecessingBatch, inputdir, Chinasupportfile, mosicmethod, isrgbn, outputdir,tarout_file
  COMPILE_OPT IDL2
  filelists=FILE_SEARCH(inputdir,'*.xml',count = num)
  PRINT,'共找到【'+STRING((num)/2, format='(i)')+'】景影像'
  PRINT,'********************************************'
  FOR i=0, num-1, 2 DO BEGIN
    msfile=filelists[i]
    pnfile=filelists[i+1]
    PRINT,'正在处理第【'+string((i/2+1), format='(i)')+'】景影像'
    end_file = PrecessingSingle(msfile, pnfile, Chinasupportfile, mosicmethod, isrgbn, outputdir)
    PRINT,'完成第【'+string((i/2+1), format='(i)')+'】景影像预处理'
    PRINT,'___________________________________________'
    print,string(end_file)
    tarout_file=end_file
  ENDFOR
  RETURN,tarout_file
END

五、使用注意事项

  1. 如果需要赋值上述代码块自行创建工程文件,注意文件名与主函数名和事件名(pretreat)保持一致,命名为pretreat.pro文件;
  2. 为保证程序运行,拼装上述代码块有前后顺序要求,建议代码块的复制先后顺序为⑤④③②①;
  3. 修订后的第二版程序支持通过国产高分支持插件,能够自动读取影像的辐射定标系数等元数据,随着时间流失,该插件可能也在不断更新,届时请自行下载最新插件或仿照其他数据格式,自行添加更新插件文件中的最新年份辐射定标系数等元数据。

六、百度云资源链接

**百度云资源链接:https://pan.baidu.com/s/1dAJaRidYGqO7-uwwEeCuVA
提取码:r6sj

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐