1. 我为什么把手工重投影换成IDL脚本GLT与RST的差异去年有段时间我一直在处理多时相影像每天的工作就是在ENVI界面里点Reproject把不同传感器的数据统一到同一网格上。白天点得烦晚上还要手动记录每次重投影用的投影参数、像元大小和输出范围。连着加了几天班之后我意识到真正该改的不是操作频率而是操作方式。于是我把整套流程搬到了IDL脚本里直接调用ENVI的GLT接口用几十行代码替代了重复劳动顺带把之前手工操作时很难发现的配准误差也暴露了出来。如果你也被投影不一致折磨过或者手头有多源影像要对齐那这篇内容应该能帮你省下不少时间。先说说核心概念GLT是Geographic Lookup Table的缩写中文经常叫地理查找表。它和传统RST多项式重投影是两条完全不同的技术路线理解清楚了后面写IDL脚本时才不会抓瞎。1.1 RST多项式重投影的近似本质RSTRotation/Scaling/Translation这种思路本质上是用一个全局多项式去拟合输入影像和输出影像之间的坐标映射关系。操作上通常是选控制点——要么手工选要么用已有几何信息——然后解算一组多项式系数最后用这套系数逐像元重采样。这套方法的好处是快、模型简单对平坦区域和覆盖范围不大的影像效果完全够用。但问题也出在全局拟合上当影像覆盖范围变大、地形起伏明显或者传感器本身的扫描几何比较复杂时一个低阶多项式很难描述那种局部畸变。就算你把多项式阶数提高也只是在控制点上拟合得更好控制点之外的地方照样可能歪。做多时相严格对比的时候这种误差会污染变化检测结果产生很多假变化。1.2 GLT逐像元查找表的工作机制GLT的思路完全不同。它不试图用一个公式去概括整幅影像的变形规律而是先把输出影像目标网格上每一个像元通过严格的投影变换反算到输入影像上应该取哪个位置的值。这个计算结果保存成一张查找表——所以叫Geographic Lookup Table。ENVI生成的GLT文件通常是两个波段一个波段存输出像元对应的输入行号另一个存输入列号。0值表示该输出像元在输入影像上没有对应位置也就是落在原始覆盖范围外。真正执行重投影的时候只需要按这张表去查位置、做插值就可以把输出影像填出来。说白了RST是算公式GLT是查字典。查字典的好处是每个像元都能得到和实际几何相符的映射关系精度高尤其适合地形复杂、传感器视角变化大的场景。1.3 什么时候必须用GLT而不是RST我的经验是以下几类情况尽量直接上GLT第一多时相影像要做逐像元对比比如NDVI时序分析拟合误差会让结果里出现假变化第二不同传感器数据要融合或叠加它们的投影、分辨率、扫描方式都不一样靠多项式硬拟合容易错位第三研究区地形起伏明显或者影像本身是宽幅、大倾角传感器获取的RST很难压住局部变形。当然GLT不是没有代价。它要先生成查找表这步计算量比RST求系数要大不少生成的文件也占磁盘空间。但对批量任务来说同一个研究区的GLT可以生成一次、反复使用所有影像共享同一张表这样反而省时而且保证了所有输出像元严丝合缝地落在同一网格上。这也是我最终选择在脚本里驱动GLT的根本原因。2. 调用ENVI接口前的环境准备初始化、文件打开、投影信息在IDL里调用ENVI接口本质上是在IDL的运行时里加载ENVI的功能库。很多从GUI转过来的用户会忽略这一步的细节上来就想写核心计算代码结果卡在环境初始化上。这里我按自己实际操作经验把准备阶段拆开讲清楚。2.1 两种初始化方式怎么选ENVI和IDL的关系比较特殊ENVI本身就是用IDL开发的所以ENVI是以功能库的形式挂在IDL运行时里的。打开IDL控制台后文件操作、投影转换、影像算法这些能力都需要先初始化ENVI才能调用。最经典的方式是直接用ENVI_INITENVI_INIT, /RESTORE_BASE_SAVE_FILES, /ENVI_DEF_FILES两个关键字分别表示恢复基础保存文件、加载默认的ENVI定义文件。如果你用的是ENVI 5.x之后的版本也可以走面向对象的方式比如e ENVI()拿到ENVI应用对象后再操作。但经典Doit系列函数在很长一段时间里都没有被移除网上沉淀下来的代码和解决方案也多所以我给工程的批量处理脚本基本还是走经典接口稳定、可迁移性强。这篇内容里的示例也以经典接口为主。2.2 打开文件与DIMS的含义初始化之后第一步是打开输入影像ENVI_OPEN_FILE, in_file, R_FIDfid IF (fid EQ -1) THEN MESSAGE, 无法打开文件: in_file ENVI_FILE_QUERY, fid, DIMSdims, NSns, NLnl, NBnb, $ DATA_TYPEdata_type, MAP_INFOmap_infoENVI_OPEN_FILE返回的是文件ID即代码里的fid。这个ID后续所有ENVI接口都会用到类似打开文件后拿到的文件句柄。ENVI_FILE_QUERY用来取文件的基本元数据包括列数ns、行数nl、波段数nb、数据类型和地图信息。这里必须重点说一下DIMS。ENVI的Doit系列函数几乎都要传DIMS参数它决定了程序处理影像的哪一块范围。最常用的全图范围写法是dims [-1, 0, ns-1, 0, nl-1]这个五元素数组的含义分别是第一个固定为-1表示按波段全部处理后面四个元素是列方向起止、行方向起止。我见过有人在批处理时图省事直接写DIMSdims但dims是之前另一景影像的结果就是输出影像要么缺一块要么在某个环节直接报错。每打开一个文件都应该用FILE_QUERY重新取一次DIMS。2.3 投影信息的读取与目标投影创建重投影最核心的输入之一是输入影像当前使用的投影。如果影像文件里已经带地图信息ENVI_FILE_QUERY返回的MAP_INFO里就包含了投影描述。我通常用ENVI_GET_PROJECTION把投影信息单独取出来proj_src ENVI_GET_PROJECTION(FIDfid)这样得到的proj_src是一个IDL结构描述当前坐标系统。目标投影一般不需要从文件读而是按需创建。以最常见的UTM投影为例proj_dst MAP_PROJ_INIT(UTM, ZONE50, DATUMWGS-84, UNITS1)UNITS1表示坐标单位为米DATUMWGS-84是基准面。不同地区、不同用途需要的投影不同有的是Albers等积投影有的是横轴墨卡托。没有把握的时候最好先找到一幅已经是你想要网格的参考影像直接读取它的投影结构比我这样手工创建要稳得多。3. 生成GLT查找表的完整步骤输出网格、边界计算、核心代码环境准备好之后就到了GLT重投影最关键的一步生成查找表。这一步的成败很大程度上取决于输出网格的设计。网格设计错了后面跑出来的影像不是全黑就是严重错位。3.1 输出网格的三个要素先想清楚不管用什么方式生成GLT你都得先确定目标网格的三要素左上角地图坐标X0、Y0、像元大小、行列数。你可以把它们理解为一张空白画布的定位和尺寸左上角坐标告诉程序画布放哪里像元大小决定分辨率行列数决定画布有多宽多高。实际项目里最省心也最常用的做法是对齐参考影像——找一幅已经位于目标网格上的影像比如某个公共底图或前一期处理好的产品直接把它左上角坐标、像元大小、行列数抄过来。这样GLT生成后输出影像和参考影像的像元能一一对应做时间序列分析时非常舒服。如果没有参考影像就得自己算。需要先获取输入影像四角的地图坐标用投影变换把它们转换到目标投影坐标系取最小外接矩形作为输出范围再根据期望的像元大小推算出行列数。这个过程容易踩半像元偏移的坑我在后面第五部分会细说。3.2 关键参数与边界计算的注意事项生成GLT时ENVI需要知道的目标投影与输出网格参数通常包括X0、Y0、PIXEL_SIZE和PROJ。X0、Y0是输出影像左上角在目标投影坐标下的位置PIXEL_SIZE是输出像元尺寸。这里有一个我很早就忽略的细节不同投影坐标系下像元大小和坐标单位的含义不同。如果目标投影坐标单位是米那么PIXEL_SIZE30就表示30米分辨率如果你不小心把输出投影设成了地理坐标系坐标单位是度那PIXEL_SIZE30就变成了30度显然荒诞。用参考影像对齐的时候我通常会直接从参考影像的地图信息里取这些参数。参考影像的MAP_INFO里包含了左上角坐标和像元尺寸直接传给GLT生成函数。要注意的是参考影像自己如果带的是斜平面投影或者特殊坐标系直接复用反而会把问题传染到新影像上所以选参考影像时要确认它的投影是规范的、与你想要的最终坐标系一致。3.3 一个可运行的IDL调用ENVI接口示例下面这段代码是我实际工程脚本的简化版本流程是打开输入影像、定义目标投影、用参考影像确定输出网格、生成GLT。PRO glt_lookup_create, in_file, ref_file, glt_file COMPILE_OPT idl2 ; 初始化ENVI ENVI_INIT, /RESTORE_BASE_SAVE_FILES, /ENVI_DEF_FILES ; 打开输入影像 ENVI_OPEN_FILE, in_file, R_FIDfid IF (fid EQ -1) THEN MESSAGE, 无法打开输入影像 ENVI_FILE_QUERY, fid, DIMSdims, NSns, NLnl ; 打开参考影像输出网格参数从它身上取 ENVI_OPEN_FILE, ref_file, R_FIDref_fid IF (ref_fid EQ -1) THEN MESSAGE, 无法打开参考影像 ENVI_FILE_QUERY, ref_fid, NSref_ns, NLref_nl, MAP_INFOref_mapinfo ; 从参考影像地图信息中取出目标网格三要素 ; 不同ENVI版本的MAP_INFO字段名可能不同用HELP, ref_mapinfo确认 x0 ref_mapinfo.x0 y0 ref_mapinfo.y0 ps ref_mapinfo.ps ; 目标投影直接用参考影像的投影最稳 proj_dst ENVI_GET_PROJECTION(FIDref_fid) ; 生成GLT查找表 ENVI_DO_GLT, FIDfid, DIMSdims, $ X0x0, Y0y0, PIXEL_SIZEps, PROJproj_dst, $ OUT_NAMEglt_file, R_FIDglt_fid PRINT, GLT已生成: , glt_file ; 清理文件句柄 ENVI_FILE_MNG, ref_fid, /REMOVE ENVI_FILE_MNG, fid, /REMOVE END这段代码的核心思路就是参考影像定网格输入影像出内容。ENVI_DO_GLT会逐像元计算目标网格上的位置对应到输入影像的哪一行哪一列结果写到glt_file里。调用完成后可以用ENVI_OPEN_FILE把GLT文件打开看一眼确认它确实是两个浮点波段。提示不同ENVI版本对相同函数的关键字命名并不完全一致。我在脚本里写的PIXEL_SIZE、X0、Y0是经典接口里的常见写法但遇到关键字不识别时在控制台输入? ENVI_DO_GLT查看你当前版本的签名比在网络上到处搜可靠得多。4. 应用GLT执行重投影插值选择、批量脚本、大影像注意GLT生成之后重投影本身反而变得很简单。ENVI提供的接口会读取查找表逐个输出像元去输入影像里取对应的值同时做插值处理。这一节把应用环节的细节和工程化注意事项讲透。4.1 ENVI_DO_GEOMETRIC调用细节与插值方式对比应用GLT做最终重投影我用的函数是ENVI_DO_GEOMETRICENVI_DO_GEOMETRIC, FIDfid, DIMSdims, GLT_FIDglt_fid, $ INTERPinterp, OUT_NAMEout_file, R_FIDout_fid关键参数就是输入文件ID和GLT文件ID。INTERP选择插值方式我习惯用0或1很少用2。三者的区别其实很直接插值方式取值特点推荐场景最近邻0直接取最近像元值速度快保持原始像元值分类图、栅格标签双线性1取周围4个像元加权平均结果平滑NDVI、连续变量三次卷积2取周围16个像元精度高但计算量大制图输出、质量要求最高这里最想提醒的是做地表覆盖分类产品重投影时千万用最近邻。我自己犯过用双线性插值处理分类结果导致边缘出现原本不存在的类别值那叫一个后悔。分类图的像元值是类别编号编号介于两个类之间没有任何物理意义。4.2 把单景流程扩展成批量脚本工程上和写Demo最大的不同是场景数量。你不可能每景影像都手动跑一遍GUI也不应该每景都重新生成GLT。同一研究区、同一目标网格的前提下GLT生成一次反复用于所有影像即可。批量脚本的骨架大概是PRO batch_glt_apply, glt_file, in_files, out_files, interp COMPILE_OPT idl2 ENVI_INIT, /RESTORE_BASE_SAVE_FILES, /ENVI_DEF_FILES ; 打开GLT文件 ENVI_OPEN_FILE, glt_file, R_FIDglt_fid IF (glt_fid EQ -1) THEN MESSAGE, GLT文件打开失败 FOR i 0, N_ELEMENTS(in_files)-1 DO BEGIN ; 打开当前影像 ENVI_OPEN_FILE, in_files[i], R_FIDfid IF (fid EQ -1) THEN CONTINUE ENVI_FILE_QUERY, fid, DIMSdims, NSns, NLnl ; 应用GLT重投影 ENVI_DO_GEOMETRIC, FIDfid, DIMSdims, GLT_FIDglt_fid, $ INTERPinterp, OUT_NAMEout_files[i], R_FIDout_fid PRINT, 完成: , in_files[i] ; 关闭本景文件防止句柄泄漏 ENVI_FILE_MNG, fid, /REMOVE IF (out_fid NE -1) THEN ENVI_FILE_MNG, out_fid, /REMOVE ENDFOR ENVI_FILE_MNG, glt_fid, /REMOVE END这里有个我自己踩过的坑批量循环里文件句柄不释放。最初我写的脚本没有在循环末尾调用ENVI_FILE_MNG结果处理到二十多景时ENVI直接拒绝再打开新文件。当时排查了很久才发现是句柄占用问题。循环内每打开一个文件及时清理就不会遇到这种莫名其妙的卡死。4.3 大影像的内存控制与中间结果验证涉及大影像时要特别注意内存。ENVI_DO_GEOMETRIC处理超大影像时如果直接全图计算内存占用可能非常夸张。更稳妥的做法是先处理一个子区验证流程正确再跑全图。你可以临时把DIMS改成一个较小的范围比如输出某一块局部dims_sub [-1, 0, 999, 0, 999] ENVI_DO_GEOMETRIC, FIDfid, DIMSdims_sub, GLT_FIDglt_fid, $ INTERPinterp, OUT_NAMEout_sub, R_FIDout_fid这个先小区后全图的习惯帮我避免了很多次跑了一整夜早上发现范围算错只能作废重来的悲剧。另外输出文件尽量用ENVI能高效读取的格式比如ENVI标准栅格格式而不是通用格式后者在大批量处理时IO开销明显更高。5. 我实测中踩过的坑和排查清单最后这部分是全文最想分享的干货。GLT重投影的代码逻辑并不复杂真正让人头疼的是生成结果看起来好像没问题实际却存在各种细微错误。我把阶段性的踩坑记录整理成清单按症状、原因、排查顺序来写。5.1 输出影像几乎全黑或者全是0值症状重投影后的影像大面积为0只有零星的像元有值。这个原因几乎可以锁定在GLT查找表没有有效覆盖输入影像。要么是输出网格的左上角坐标和像元尺寸与目标投影根本不匹配要么是PROJ参数给错了。比如输入影像本来在UTM 50N目标投影却设置成了UTM 49N跨了带以后坐标范围对不上查找表里自然全是0。排查时不要先去改重投影代码先打开GLT文件看看它的像元统计值。正常来说GLT的两个波段最大值应该约等于输入影像的行数和列数。如果最大值是0或者超过了输入尺寸那问题一定出在生成GLT那一步从头检查投影和网格参数。5.2 影像内部出现条纹或锯齿状错位症状重投影结果边缘黑边尚可接受但影像内部有规律的错位条纹。内部条纹一般不是投影定义错而是两类问题。第一DIMS不匹配。生成GLT时用的是某一景影像的DIMS应用GLT时却用了另一景的DIMS导致查找表与实际数据空间错位。第二插值方式和数据性质不匹配。连续变量用三次卷积或双线性都还好但如果你处理的是有NaN等特殊值的影像插值会把特殊值扩散开看起来就像条纹。处理方案是给GLT生成和应用环节统一套同一份DIMS来源特殊值在重投影前先掩膜输出后再把非有效区域重新标记。5.3 差半像元地图坐标中心点与角点的迷惑这个坑特别隐蔽。ENVI里MAP_INFO记录的坐标到底是像元中心还是像元左上角在不同文件、不同处理流程里含义不完全一样。当我用参考影像的MAP_INFO去驱动GLT时如果参考影像的网格本身和另一个数据存在半像元平移所有输出结果都会统一偏移肉眼可能看不出但和矢量边界套合时就会露出破绽。判断方法很土但管用在输入影像和参考影像上找同一个稳定的地物点比如道路交叉口比较它在两个结果里的像素坐标。如果偏差稳定在一个像元左右那就检查网格原点的定义按半像元修正X0、Y0。这个土办法比任何理论分析都快。5.4 函数签名版本差异引起的接口报错IDL调用ENVI接口的代码网上能找到的很多都是多年前的版本。ENVI升级后部分关键字被改过名比如我早期就遇到过PIXEL_SIZE在某些版本中变成了O_PIXEL_SIZE。如果你从网上抄了一段代码运行时报Unknown keyword之类的错误不要怀疑自己理解错了先查当前环境的签名。查签名的方式很简单在IDL控制台输入IDL ? ENVI_DO_GLT这会打开当前版本对该函数参数的详细说明。另外HELP, /ROUTINES, ENVI_DO_GLT也可以确认函数是否在当前环境被加载。我踩过几次版本坑之后养成了代码从网上拿来先查签名再改参数的习惯基本没有再被这类问题卡住过。5.5 一个十秒钟的自查习惯能挡掉一大半错误最后分享一个我每次跑完GLT必做的小动作打开生成的GLT文件点开像元浏览随便找几个不同的区域看看波段值。只要看到行列号的最大值接近输入影像的尺寸、而且分布有层次没有大面积0值我就知道GLT没生成废。这个习惯只需要十秒钟却已经帮我挡掉了无数次因为网格参数笔误导致的返工。整个过程走下来你会发现IDL调用ENVI做GLT重投影并没有多玄妙难的是把参数的含义真正吃透以及培养一种先查表、再重投影、最后验结果的工程习惯。数据越多、时间越紧的时候这种脚本化、流程化的处理方式越能体现价值。