1. 我为什么把手工重投影换成IDL脚本：GLT与RST的差异
去年有段时间我一直在处理多时相影像，每天的工作就是在ENVI界面里点Reproject，把不同传感器的数据统一到同一网格上。白天点得烦，晚上还要手动记录每次重投影用的投影参数、像元大小和输出范围。连着加了几天班之后我意识到，真正该改的不是操作频率，而是操作方式。于是我把整套流程搬到了IDL脚本里，直接调用ENVI的GLT接口，用几十行代码替代了重复劳动，顺带把之前手工操作时很难发现的配准误差也暴露了出来。
如果你也被投影不一致折磨过，或者手头有多源影像要对齐，那这篇内容应该能帮你省下不少时间。先说说核心概念：GLT是Geographic Lookup Table的缩写，中文经常叫地理查找表。它和传统RST多项式重投影是两条完全不同的技术路线，理解清楚了，后面写IDL脚本时才不会抓瞎。
1.1 RST多项式重投影的近似本质
RST（Rotation/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_INIT：
ENVI_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_FID=fid IF (fid EQ -1) THEN MESSAGE, '无法打开文件: ' + in_file ENVI_FILE_QUERY, fid, DIMS=dims, NS=ns, NL=nl, NB=nb, $ DATA_TYPE=data_type, MAP_INFO=map_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，表示按波段全部处理；后面四个元素是列方向起止、行方向起止。我见过有人在批处理时图省事，直接写DIMS=dims但dims是之前另一景影像的，结果就是输出影像要么缺一块，要么在某个环节直接报错。每打开一个文件，都应该用FILE_QUERY重新取一次DIMS。
2.3 投影信息的读取与目标投影创建
重投影最核心的输入之一，是输入影像当前使用的投影。如果影像文件里已经带地图信息，ENVI_FILE_QUERY返回的MAP_INFO里就包含了投影描述。我通常用ENVI_GET_PROJECTION把投影信息单独取出来：
proj_src = ENVI_GET_PROJECTION(FID=fid)这样得到的proj_src是一个IDL结构，描述当前坐标系统。目标投影一般不需要从文件读，而是按需创建。以最常见的UTM投影为例：
proj_dst = MAP_PROJ_INIT('UTM', ZONE=50, DATUM='WGS-84', UNITS=1)UNITS=1表示坐标单位为米，DATUM='WGS-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_SIZE=30就表示30米分辨率；如果你不小心把输出投影设成了地理坐标系，坐标单位是度，那PIXEL_SIZE=30就变成了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_FID=fid IF (fid EQ -1) THEN MESSAGE, '无法打开输入影像' ENVI_FILE_QUERY, fid, DIMS=dims, NS=ns, NL=nl ; 打开参考影像，输出网格参数从它身上取 ENVI_OPEN_FILE, ref_file, R_FID=ref_fid IF (ref_fid EQ -1) THEN MESSAGE, '无法打开参考影像' ENVI_FILE_QUERY, ref_fid, NS=ref_ns, NL=ref_nl, MAP_INFO=ref_mapinfo ; 从参考影像地图信息中取出目标网格三要素 ; 不同ENVI版本的MAP_INFO字段名可能不同，用HELP, ref_mapinfo确认 x0 = ref_mapinfo.x0 y0 = ref_mapinfo.y0 ps = ref_mapinfo.ps ; 目标投影：直接用参考影像的投影，最稳 proj_dst = ENVI_GET_PROJECTION(FID=ref_fid) ; 生成GLT查找表 ENVI_DO_GLT, FID=fid, DIMS=dims, $ X0=x0, Y0=y0, PIXEL_SIZE=ps, PROJ=proj_dst, $ OUT_NAME=glt_file, R_FID=glt_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_GEOMETRIC：
ENVI_DO_GEOMETRIC, FID=fid, DIMS=dims, GLT_FID=glt_fid, $ INTERP=interp, OUT_NAME=out_file, R_FID=out_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_FID=glt_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_FID=fid IF (fid EQ -1) THEN CONTINUE ENVI_FILE_QUERY, fid, DIMS=dims, NS=ns, NL=nl ; 应用GLT重投影 ENVI_DO_GEOMETRIC, FID=fid, DIMS=dims, GLT_FID=glt_fid, $ INTERP=interp, OUT_NAME=out_files[i], R_FID=out_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, FID=fid, DIMS=dims_sub, GLT_FID=glt_fid, $ INTERP=interp, OUT_NAME=out_sub, R_FID=out_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重投影并没有多玄妙，难的是把参数的含义真正吃透，以及培养一种"先查表、再重投影、最后验结果"的工程习惯。数据越多、时间越紧的时候，这种脚本化、流程化的处理方式越能体现价值。
