基于IDL的高分二号影像批量预处理程序
基于IDL的高分二号遥感影像批量预处理程序(附件含源码和sav文件),预处理包括解压、几何校正、辐射定标、大气校正、自动校正、融合。
基于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:是否包含近红外波段(可选值有1和0);
;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:是否包含近红外波段(可选值有1和0);
;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
五、使用注意事项
- 如果需要赋值上述代码块自行创建工程文件,注意文件名与主函数名和事件名(pretreat)保持一致,命名为pretreat.pro文件;
- 为保证程序运行,拼装上述代码块有前后顺序要求,建议代码块的复制先后顺序为⑤④③②①;
- 修订后的第二版程序支持通过国产高分支持插件,能够自动读取影像的辐射定标系数等元数据,随着时间流失,该插件可能也在不断更新,届时请自行下载最新插件或仿照其他数据格式,自行添加更新插件文件中的最新年份辐射定标系数等元数据。
六、百度云资源链接
**百度云资源链接:https://pan.baidu.com/s/1dAJaRidYGqO7-uwwEeCuVA
提取码:r6sj
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)