高分二号
下面写的太LOW了，有好多值得修改的地方，修改更新后博客地址：
https://blog.csdn.net/desertsTsung/article/details/84679969
今天去不了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
; [email protected]
PRO 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.Close
END
文章最后发布于: 2018-11-22 00:44:58
相关阅读
文章对各大电商平台的收货地址进行了系列的对比分析，尝试归纳出一套体验较好、操作较快，认知明确的收货地址流程。希望能够本文对你
turbosun是甄嬛娘娘孙俪的微博ID，近日自曝开始做起了老板在网上卖家里的闲置物品，虽然吐槽太辛苦要改“邪”归正好好拍戏，但却揭开
实现自动化安装操作系统我们仍需要插入光盘来引导，现在很多服务器已经没有光驱，那么此时我们就无法用光盘引导，如果要实现光盘引导安
A5创业网(公众号：iadmin5)11月15日报道，继113个10000元现金抽奖后，王思聪又开启了抽奖第二波，转发微博送67套纪念版手机壳，为了防止型
淘宝链接