本文共 4972 字，大约阅读时间需要 16 分钟。

下面写的太LOW了，有好多值得修改的地方，修改更新后博客地址：

今天去不了Earth Engine，所以就写了一点IDL的代码。

主要包括对高分二号数据的正射校正，（裁剪），辐射定标，快速大气校正，影像自动配准和影像融合。

如果用ENVI来做如上预处理的话会产生一些临时文件需要自己删除，而且每个过程都需要自己亲自操作，处理时间过于碎片化。基于以上原因写了文末的代码，可以自己输入高分二号的文件夹名来一步直达最后的融合结果，减少了中间的人为干预，也不需要自己删中间文件，可以说是很大的优化了处理时间。结果图：

说明：代码中用到了很多ENVITask，每个Task引入IDL的版本都不一样，所以最低版本要求为ENVI5.2.1。16行dmfile为正射校正的DEM文件，代码中默认为ENVI自带的全球DEM文件，如果对精度有要求，可以自行修改。至于时间，1/64的影像耗时在半小时以下，如果未裁剪，保守估计所有步骤跑下来应该在一个小时以上。  实测中一景影像耗时为260分钟左右，测试电脑配置为i5-6200U，4GB内存。如下源码为处理整个一景。

;+; Description:;    Pre-Treatment of GaoFen-2 Imagery;; Author: JAYTON TSUNG Nov 21 ,2018;         CHENGDU UNIVERSITY OF INFORMATION TECHNOLOGY;         DESERTSTSUNG@QQ.COMPRO PreTreat  COMPILE_OPT IDL2  e = ENVI(/headless)  t = SYSTIME(1)    indir = 'H:\GF2_PMS2_E115.6_N38.4_20170628_L1A0002448743'  dmfile = e.Root_Dir + 'data\GMTED2010.jp2'  pos = STRSPLIT(indir, '\')  outdir = STRMID(indir, pos[0], pos[-1])    CD,indir  filelist = FILE_SEARCH('*.tiff')  msfile = filelist[0] & pnfile=filelist[1]    ;1_RPCOrthorectification  PRINT,'Part1 of 5 Has Already Started, Please Wait'  ort = ENVITask('RPCOrthorectification') ;ENVI5.1 Introduced  ms = e.Openraster(msfile) & dm = e.Openraster(dmfile)  ort.dem_raster = dm  ort.input_raster = ms  ort.Output_Raster_URI = e.GetTemporaryFilename()  ort.EXECUTE    ortpn = ENVITask('RPCOrthorectification')  pn = e.Openraster(pnfile) & dm = e.Openraster(dmfile)  ortpn.dem_raster = dm  ortpn.input_raster = pn  ortpn.Output_Raster_URI = e.GetTemporaryFilename()  ortpn.EXECUTE  PRINT,'Part1 of 5 Finished, Takes '$    + STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes'    ;2_ApplyGainOffset  PRINT,'Part2 of 5 Has Already Started, Please Wait'  IF STRPOS(indir, 'PMS1') NE -1 THEN BEGIN ;;;;PMS1;;;;    IF STRPOS(indir, '2015') NE -1 THEN BEGIN      msgain = [0.1457, 0.1604, 0.155, 0.1731]      pngain = [0.1538]      msoffs = [0, 0, 0, 0]      pnoffs = [0]    ENDIF ELSE IF STRPOS(indir, '2016') NE -1 THEN BEGIN      msgain = [0.1322, 0.155, 0.1477, 0.1613]      pngain = [0.1501]      msoffs = [0, 0, 0, 0]      pnoffs = [0]    ENDIF ELSE IF STRPOS(indir, '2017') NE -1 THEN BEGIN      msgain = [0.1193, 0.153, 0.1424, 0.1569]      pngain = [0.1503]      msoffs = [0, 0, 0, 0]      pnoffs = [0]    ENDIF  ENDIF ELSE BEGIN ;;;;PMS2;;;;    IF STRPOS(indir, '2015') NE -1 THEN BEGIN      msgain = [0.1761, 0.1843, 0.1677, 0.183]      pngain = [0.1538]      msoffs = [0, 0, 0, 0]      pnoffs = [0]    ENDIF ELSE IF STRPOS(indir, '2016') NE -1 THEN BEGIN      msgain = [0.1762, 0.1856, 0.1754, 0.1980]      pngain = [0.1863]      msoffs = [0, 0, 0, 0]      pnoffs = [0]    ENDIF ELSE IF STRPOS(indir, '2017') NE -1 THEN BEGIN      msgain = [0.1434, 0.1595, 0.1511, 0.1685]      pngain = [0.1679]      msoffs = [0, 0, 0, 0]      pnoffs = [0]    ENDIF  ENDELSE    rad = ENVITask('ApplyGainOffset') ;ENVI5.2.1 Introduced  ms = ort.Output_Raster  rad.Input_Raster = ms  rad.Gain = msgain & rad.offset = msoffs  rad.Output_Raster_URI = e.GetTemporaryFilename()  rad.EXECUTE    radpn = ENVITask('ApplyGainOffset')  pn = ortpn.Output_Raster  radpn.Input_Raster = pn  radpn.Gain = pngain & radpn.offset = pnoffs  radpn.Output_Raster_URI = e.GetTemporaryFilename()  radpn.EXECUTE  PRINT,'Part2 of 5 Finished, Takes '$    + STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes'    ;3_QUAC  PRINT,'Part3 of 5 Has Already Started, Please Wait'  ms = rad.Output_Raster  mtd = ms.Metadata ;ENVI5 Introduced  mtd.AddItem, 'Wavelength', [492,555,665,822]  mtd.AddItem, 'Wavelength Units', 'Nanometers'    qac = ENVITask('QUAC') ;ENVI5.1 Introduced  ms = rad.Output_Raster  qac.Input_Raster = ms  qac.Output_Raster_URI = e.GetTemporaryFilename()  qac.EXECUTE  PRINT,'Part3 of 5 Finished, Takes '$    + STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes'    ;4_ImageToImageRegistration  PRINT,'Part4 of 5 Has Already Started, Please Wait'  tps = ENVITask('GenerateTiePointsByCrossCorrelation') ;ENVI5.2.1 Introduced  ms = qac.Output_Raster  pn = radpn.Output_Raster  tps.Input_Raster1 = pn  tps.Input_Raster2 = ms  tps.EXECUTE    flt = ENVITask('FilterTiePointsByGlobalTransform') ;ENVI5.2.1 Introduced  TiePoints = tps.Output_Tiepoints  flt.Input_Tiepoints = TiePoints  flt.EXECUTE    rgs = ENVITask('ImageToImageRegistration') ;ENVI5.2.1 Introduced  TiePoints2 = flt.Output_Tiepoints  rgs.Input_Tiepoints = TiePoints2  rgs.Warping = 'Triangulation'  rgs.Output_Raster_URI = e.GetTemporaryFilename()  rgs.EXECUTE  PRINT,'Part4 of 5 Finished, Takes '$    + STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes'    ;5_GramSchmidtPanSharpening  PRINT,'Part5 of 5 Has Already Started, Please Wait'  gsf = ENVITask('GramSchmidtPanSharpening') ;ENVI5.2 Introduced  ms = rgs.Output_Raster  pn = radpn.Output_Raster  gsf.Input_Low_Resolution_Raster = ms  gsf.Input_High_Resolution_Raster = pn  gsf.Output_Raster_URI = outdir + 'GS_Fusion'  gsf.EXECUTE  PRINT,'Part5 of 5 Finished, Takes '$    + STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes'    e.CloseEND