IDL下高分二号完整预处理代码

下面写的太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
;         DESERTSTSUNG@QQ.COMPRO PreTreatCOMPILE_OPT IDL2e = 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,indirfilelist = FILE_SEARCH('*.tiff')msfile = filelist[0] & pnfile=filelist[1];1_RPCOrthorectificationPRINT,'Part1 of 5 Has Already Started, Please Wait'ort = ENVITask('RPCOrthorectification') ;ENVI5.1 Introducedms = e.Openraster(msfile) & dm = e.Openraster(dmfile)ort.dem_raster = dmort.input_raster = msort.Output_Raster_URI = e.GetTemporaryFilename()ort.EXECUTEortpn = ENVITask('RPCOrthorectification')pn = e.Openraster(pnfile) & dm = e.Openraster(dmfile)ortpn.dem_raster = dmortpn.input_raster = pnortpn.Output_Raster_URI = e.GetTemporaryFilename()ortpn.EXECUTEPRINT,'Part1 of 5 Finished, Takes '$+ STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes';2_ApplyGainOffsetPRINT,'Part2 of 5 Has Already Started, Please Wait'IF STRPOS(indir, 'PMS1') NE -1 THEN BEGIN ;;;;PMS1;;;;IF STRPOS(indir, '2015') NE -1 THEN BEGINmsgain = [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 BEGINmsgain = [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 BEGINmsgain = [0.1193, 0.153, 0.1424, 0.1569]pngain = [0.1503]msoffs = [0, 0, 0, 0]pnoffs = [0]ENDIFENDIF ELSE BEGIN ;;;;PMS2;;;;IF STRPOS(indir, '2015') NE -1 THEN BEGINmsgain = [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 BEGINmsgain = [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 BEGINmsgain = [0.1434, 0.1595, 0.1511, 0.1685]pngain = [0.1679]msoffs = [0, 0, 0, 0]pnoffs = [0]ENDIFENDELSErad = ENVITask('ApplyGainOffset') ;ENVI5.2.1 Introducedms = ort.Output_Rasterrad.Input_Raster = msrad.Gain = msgain & rad.offset = msoffsrad.Output_Raster_URI = e.GetTemporaryFilename()rad.EXECUTEradpn = ENVITask('ApplyGainOffset')pn = ortpn.Output_Rasterradpn.Input_Raster = pnradpn.Gain = pngain & radpn.offset = pnoffsradpn.Output_Raster_URI = e.GetTemporaryFilename()radpn.EXECUTEPRINT,'Part2 of 5 Finished, Takes '$+ STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes';3_QUACPRINT,'Part3 of 5 Has Already Started, Please Wait'ms = rad.Output_Rastermtd = ms.Metadata ;ENVI5 Introducedmtd.AddItem, 'Wavelength', [492,555,665,822]mtd.AddItem, 'Wavelength Units', 'Nanometers'qac = ENVITask('QUAC') ;ENVI5.1 Introducedms = rad.Output_Rasterqac.Input_Raster = msqac.Output_Raster_URI = e.GetTemporaryFilename()qac.EXECUTEPRINT,'Part3 of 5 Finished, Takes '$+ STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes';4_ImageToImageRegistrationPRINT,'Part4 of 5 Has Already Started, Please Wait'tps = ENVITask('GenerateTiePointsByCrossCorrelation') ;ENVI5.2.1 Introducedms = qac.Output_Rasterpn = radpn.Output_Rastertps.Input_Raster1 = pntps.Input_Raster2 = mstps.EXECUTEflt = ENVITask('FilterTiePointsByGlobalTransform') ;ENVI5.2.1 IntroducedTiePoints = tps.Output_Tiepointsflt.Input_Tiepoints = TiePointsflt.EXECUTErgs = ENVITask('ImageToImageRegistration') ;ENVI5.2.1 IntroducedTiePoints2 = flt.Output_Tiepointsrgs.Input_Tiepoints = TiePoints2rgs.Warping = 'Triangulation'rgs.Output_Raster_URI = e.GetTemporaryFilename()rgs.EXECUTEPRINT,'Part4 of 5 Finished, Takes '$+ STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes';5_GramSchmidtPanSharpeningPRINT,'Part5 of 5 Has Already Started, Please Wait'gsf = ENVITask('GramSchmidtPanSharpening') ;ENVI5.2 Introducedms = rgs.Output_Rasterpn = radpn.Output_Rastergsf.Input_Low_Resolution_Raster = msgsf.Input_High_Resolution_Raster = pngsf.Output_Raster_URI = outdir + 'GS_Fusion'gsf.EXECUTEPRINT,'Part5 of 5 Finished, Takes '$+ STRING((SYSTIME(1)-t)/60, format='(f7.2)') + ' Minutes'e.CloseEND

IDL下高分二号完整预处理代码相关推荐

ENVI5.3下高分二号（GF2）数据预处理
以一景2015年1月23日获取的GF2-PMS1数据为例介绍在ENVI5.3下GF2数据预处理的详细操作步骤.GF2数据预处理基本流程如下: 图:GF2数据预处理流程说明:1. 针对不同的应用,有不 ...
基于IDL的高分二号影像批量预处理程序
基于IDL的高分二号影像批量预处理程序(第二版) 上一版本插件预处理效果不太理想(特别是融合后分辨率不符合预期),上一版资源不再提供下载,但处理思路可以借鉴.目前在本机40GB RAM工作站上预处理一 ...
第059篇：高分二号遥感影像预处理流程（ENVI5.3.1平台+ENVI App Store中最新的中国国产卫星支持工具）
今天被袁老的新闻刷屏,湖南衡水县水稻基地传出好消息: 袁隆平团队第三代杂交水稻测产,测得晚稻平均亩产为911.7 ...
遥感图像语义分割——从原始图像开始制作自己的数据集(以高分二号为例)
遥感图像语义分割--从原始图像开始制作自己的数据集(以高分二号为例) 文章目录遥感图像语义分割--从原始图像开始制作自己的数据集(以高分二号为例) 1.遥感影像获取 2.遥感数据预处理(影像融合) ...
高分二号多光谱GEOTIFF影像批量转换为JPG图片
最近在做语义分割方面的研究,想对比一下高分二号4个波段的模型输入和JPG输入的结果有什么差别.使用了一些图片转换软件发现转换后的jpg图片一片漆黑,应该是geitif格式的特殊性吧!所有在网上找了一些 ...
ENVI高分一号卫星预处理：以GF-1/WFV多光谱数据为例（含定标系数）
目录 GF1/WFV卫星简介与预处理流程: 1.辐射定标: 2.大气校正: 3.正射/几何校正: 4.配准:具体可参考链接 5.拼接: GF1/WFV卫星简介与预处理流程: GF1/WFV卫星数据参数 ...
基于python实现高分二号遥感影像水体提取与水质反演（黑臭水体与水体富营养化）
高分二号遥感影像水体提取与水质反演水体提取函数--NDWI 基于几何约束提取河流生成shp,方便后续裁剪水体水质反演最终结果水体提取函数--NDWI 水体提取函数water.py impor ...
使用ERDAS对国产卫星影像进行控制点正射校正—以高分二号卫星（GF2）为例
** 使用ERDAS对国产卫星影像进行控制点正射校正-以高分二号卫星(GF2)为例 ** 一.研究区影像: 1.1.高分二号原始影像影像为山西介休市的一景获取时间为2019.01.22的GF2影像: ...
遥感技术及高分遥感影像在地震中的应用及高分二号获取
长期以来,地震预报监测.灾害调查.灾情信息获取主要依靠实地勘测手段,其获取的数据精度和置信度虽然较高,但存在工作量大.效率低.费用高和信息不直观等缺点.遥感技术手段可在一定程度上克服传统实地勘测手段的 ...

IDL下高分二号完整预处理代码

IDL下高分二号完整预处理代码相关推荐

最新文章

热门文章