做遥感数据处理这个方向MODIS 数据几乎是绕不开的一类源数据而只要碰过 MODIS正弦曲线投影Sinusoidal Projection和 ENVI 下的坐标投影定义这两个词就一定在某个深夜折磨过你。我自己最开始接手一批 MODIS 反射率产品时打开文件一看图像是“躺”在那儿的但坐标栏里空空如也量出来的数值一会儿是几千、一会儿是几百万完全不知道自己在看地球的哪个角落。折腾了小半天才明白MODIS 的大部分格网产品天生带着正弦曲线投影信息但这份信息藏在 HDF 的元数据结构里ENVI 有时能读到有时读不到读不到的时候就只能手工把投影定义补回去。而“补回去”这三个字背后涉及投影类型、基准球体半径、中央经线、像素大小、左上角坐标一整套参数的核对。这篇内容就围绕这条线展开正弦曲线投影到底是个什么几何构造、MODIS 为什么选中它、ENVI 里定义坐标投影的具体操作路径、参数怎么算、重投影和 GLT 校正怎么选以及我踩过的那些坑。适合已经有基础 GIS 概念、正在处理 MODIS 格网产品的同学也适合刚接触遥感影像坐标系统、想把这套逻辑一次搞明白的新手。1. 正弦曲线投影到底在解决什么问题1.1 一张“压扁的橘子皮”正弦投影的几何直觉地理坐标是球面的屏幕和数组是平面的把一个球面摊平这件事数学上注定要付出代价——要么面积变形要么形状变形要么距离变形三者不可兼得。正弦曲线投影的思路很朴素想象地球是一个标准球体沿着赤道切一刀把所有纬线按各自的周长“收拢”到一条中央经线上然后再展开。赤道上纬线最长展开后的横向跨度最大越往高纬度走纬线圈越小横向跨度按纬度余弦值收缩到极点收缩成一个点。落到公式上就是 x R·(λ − λ₀)·cosφy R·φ。R 是球体半径λ 是经度λ₀ 是中央经线φ 是纬度。你能从公式里直接读出两件事一是横向坐标被 cosφ 加权了这意味着不同纬度带上的单位长度在图上代表的实际距离不同但图上每一块面积对应的地表面积是相等的二是纵向坐标只跟纬度线性相关纬度和 y 值是一一对应的。这种“等积”特性是正弦投影最核心的价值。用生活里的说法它就像把橘子皮剥下来以后不是硬拉平而是顺着经线方向小心地压扁让每块皮的面积保持原样但形状被扭曲了。低纬度地区看起来还比较正常高纬度地区的形状会被横向拉得很夸张这也是为什么你在正弦投影的全球图上看到某些高纬度陆块会觉得“胖得离谱”。1.2 为什么 MODIS 这类产品偏爱它MODIS 的产品体系里格网产品Gridded L3占了很大比重包括反射率、植被指数、地表温度、土地覆盖等。选择正弦投影我理解是几个需求叠加的结果。第一是等积性。这类产品大量用于区域乃至全球尺度的定量分析比如统计某个纬度带的植被指数均值、比较不同区域的能量收支。如果投影本身面积不等跨纬度比较就会出现系统性偏差后续修正非常麻烦。正弦投影天然等积省掉了这一层顾虑。第二是全球无缝拼接。MODIS 的正弦投影格网把全球切成 36 列 × 18 行的瓦片tile每块覆盖经度 10°、纬度 10°标准尺寸是 1200×1200 像素1 公里产品。相邻瓦片之间在网格坐标上严丝合缝拼接时不需要做任何重叠裁剪或边缘羽化直接按行列号摆好就行。这一点在实际业务里价值巨大——你处理全球数据时不用为每块瓦片计算复杂的地理范围。第三是计算简单。正弦投影只依赖球体半径和中央经线两个参数没有标准纬线、没有假东假北、没有复杂的椭球参数读取和写入的元数据量很小对于动辄几十 GB 的全球产品来说这种轻量化的定义方式很友好。当然代价也明显高纬度形状失真严重不适合做局部尺度的制图或导航跨大范围的图上距离没有实际意义。所以正弦投影更像是“内部存储格式”最终交付或做局地分析时一般还是要转到地理坐标或适合当地的投影上。1.3 三个必须记住的投影参数不管后面用什么工具MODIS 正弦投影的定义里绕不开三个数这三个数记错了后面全白搭。参数取值说明球体半径 R6371007.181 m不是 WGS84 椭球是专门定义的球体中央经线 λ₀0°全球格网统一以本初子午线为中心全球横向范围−20015109.354 m ~ 20015109.354 m赤道处半周长球体半径这个值是关键中的关键。很多人下意识填 6378137WGS84 椭球的赤道半径结果就是整幅图偏出去几十米到上百米做小比例尺看不出来一旦和矢量边界叠加或者做局地裁剪偏移就非常刺眼。6371007.181 这个数字是 MODIS 正弦格网专门定义的平均球体半径必须原样使用。中央经线为 0意味着 x 0 的位置落在本初子午线上全球范围的 x 值从负两千万到正两千万。全球纵向范围则是从赤道往北、往南各约一千万米对应 y 从 10007554.677 到 −10007554.677 左右。2. MODIS 数据里的投影信息藏在哪2.1 HDF 元数据结构里的“自白书”MODIS 的 HDF-EOS 文件结构里有一块专门存放格网结构定义的区域通常叫 StructMetadata。这里面会明确写出投影类型、投影参数、格网的行列数、左上角坐标、右下角坐标等。如果你用 HDF 查看工具打开这个结构能读到类似ProjectionGCTP_SINUSOIDAL这样的字样后面跟着一组数字参数第一个数就是球体半径 6371007.181其余位置补零再往下是UpperLeftPointMtrs、LowerRightMtrs、XDim、YDim这类格网几何信息。这套元数据设计得很完整理论上只要软件正确解析投影定义就不需要人工干预。问题出在“正确解析”这四个字上。不同软件、不同版本、不同产品类型对 HDF-EOS 元数据的兼容程度不一样有的产品 StructMetadata 完整有的简化过有的把关键字段放在了别的属性里。我一般的习惯是拿到一批陌生的 MODIS 数据先不急着在 ENVI 里打开看影像而是先看一眼元数据结构。ENVI 5.x 里可以用 File 菜单下的元数据查看功能也可以在打开 HDF 时选择只读元信息。确认里面确实有正弦投影和格网参数再决定后续是不是需要手工定义。这一步花不了两分钟能省掉后面半小时的排查。2.2 为什么 ENVI 打开后经常读到“没有投影”实际遇到 ENVI 读不到投影原因大致能归成几类。一类是产品本身就没有写完整投影信息比如某些经过二次处理的衍生数据、别人转存过的版本StructMetadata 被裁剪或者丢失一类是 HDF 里存在多个数据集子数据集、科学数据集ENVI 打开时默认选中了某个图层但投影信息绑定在格网结构层而不是图层本身导致属性没带过来还有一类是文件被多次转存比如从 HDF 转成 TIFF 再转回来元数据在转换中被截断。还有一种容易被忽略的情况数据其实带着投影但 ENVI 把坐标单位识别成了“米”却把投影类型识别成“未定义”界面上表现为 Map Info 里有一堆数字但没有投影名称。这种情况下图像的位置其实是对的只是缺少投影标识重投影时无法判断源坐标系。判断方法很简单打开头文件看 Map Info 那一行。如果一整行是空的或者只有{0,0,0,0,0,0}之类的占位值那就是彻底没有如果有 7 个以上数字且第一个数字不是 0说明可能带着某种定义只是名称没解析出来这时候可以对照投影代码表去核对。提示不要看到图像显示正常就默认为投影正确。ENVI 在完全没有投影信息时会按像素行列号当坐标显示看起来一样“正常”但那是假坐标。2.3 定义前必须核对的五个量手工补投影之前先准备好这张核对表五个量确认无误再动手否则改来改去特别浪费时间。核对项怎么确认常见错误投影类型元数据中 Projection 字段误选 Geographic Lat/Lon球体半径元数据 ProjParams 第一个值填成 WGS84 椭球半径像素大小按产品分辨率反推1 公里产品填成 250 米值左上角坐标由瓦片行列号推算少算或多算一块的偏移数据行列数打开后看波段尺寸与像素大小不匹配导致缩放错误像素大小这一项特别容易出问题。MODIS 的正弦格网像素值不是整数250 米产品是 231.656358 米500 米产品是 463.312717 米1 公里产品是 926.625433 米。这些非整数值看着别扭但它们是全球 36×18 瓦片刚好铺满的必然结果。如果你图省事填了 250、500、1000图像会轻微缩放瓦片之间拼接时会出现细缝肉眼在小范围看不清做镶嵌的时候就露馅了。3. ENVI 下定义坐标投影的完整流程3.1 从打开文件到看清 HeaderENVI 里打开 MODIS HDF 有两条路一条是从文件菜单里选对应的数据格式入口另一条是直接拖拽文件到视图里。两种方式底层逻辑一样都会去解析 HDF 结构区别是拖拽方式有时会跳过某些元数据读取步骤。打开之后第一件事不是看影像而是打开头文件编辑界面。找到 Map Info 和 Projection Info 这两项看图 3.1 描述的状态——如果 Map Info 里已经有一串坐标值先别急着高兴把数值对照一下x 应该在负两千万到正两千万之间y 在负一千万到正一千万之间。如果填的是几百、几千这种数那多半是按行列号占位说明投影彻底没读到。接着看 Projection Info。ENVI 不同版本这里的格式略有差异常见形式是以投影代码开头的一组参数。正弦投影的代码是 16后面跟球体半径和中央经线等。如果这里显示的是“Unknown”或者干脆是空数组那就要手工补。还有一种情况是打开了 HDF 里的多个图层后做图层叠加Layer Stacking操作叠加完成的新文件投影信息丢失。这是很常见的坑因为叠加过程只保留像素数组不自动继承源投影。处理办法是叠加时手动指定参考文件或者叠加完重新定义一次投影。3.2 手工写入 Map Info 与 Projection Info在头文件编辑界面里找到 Map Info 那一栏点进去可以手动填。ENVI 的 Map Info 数组结构大致是投影名称或代码、基准面代码、左上角像素的左上角 x 坐标、左上角像素的左上角 y 坐标、x 方向像素大小、y 方向像素大小后面再跟若干投影参数占位。对正弦投影来说一组典型的手工填入值是这样map info {Sinusoidal, 1, 1, -20015109.354, 10007554.677, 926.6254331, 926.6254331, 0, 0, 0, 0, 0, 0, 0, 0}拆开看第一项是投影类型第二、三项是基准面和单位相关的标识第四、五项是左上角坐标这里对应全球格网最左上角那块瓦片列号 0、行号 0的角点第六、七项是像素大小这里填的是 1 公里产品的值后面一串 0 是正弦投影用不到的标准纬线、假东假北等参数。同时在 Projection Info 里补上projection info {16, 6371007.181, 0, 0}第一项 16 是正弦投影的代码第二项是球体半径第三项是中央经线。不同 ENVI 版本对这个数组的长度要求不完全一致如果填写后报错可以先只填前两项试一下再逐步补全。我建议每次改完都立刻保存然后关掉文件重新打开看一眼确认数值没有被软件自动重置。注意手工填 Map Info 时坐标对应的是“左上角像素的外角”不是像素中心。这两个概念差半个像素在 1 公里产品上就是 463 米落到地面相当明显。3.3 用 Georeference MODIS 工具半自动处理如果 ENVI 版本里带有专门的地理定位工具可以用它来做半自动处理。这条路径的好处是不用自己算瓦片坐标工具会根据元数据或文件名里的行列号自动推算。流程大致是在工具箱里找到对应的地理定位入口选择输入文件工具会读取格网结构弹出参数确认界面核对投影类型、球体半径、像素大小确认无误后执行。执行完之后会生成一个新的输出文件投影信息完整且通常已经完成了一次正弦投影下的网格重排。这时候再打开头文件Map Info 里应该是一组完整坐标用光标查询能看到合理的 x/y 值。这条路径的局限在于版本兼容性。老版本工具箱里有独立入口新版本可能集成到了别的菜单下或者需要额外扩展。如果你的环境里找不到不要死磕回到手工定义那条路五分钟也能搞定而且过程透明出了问题好排查。3.4 像素大小与左上角坐标的推算实例这是最值得花时间弄明白的一块。全球格网被切成 36 列 18 行每块瓦片覆盖 10°×10°瓦片的左上角坐标可以直接由行列号线性推算。整幅全球图的宽度是 2 × 20015109.354 40030218.708 米切成 36 块每块宽度 1111950.520 米。这个宽度对三种分辨率都成立因为不同分辨率只是每块内部的像素数不同250 米产品每块 4800 像素500 米产品每块 2400 像素1 公里产品每块 1200 像素乘回去都等于 1111950.520。产品分辨率单瓦片像素数像素大小米瓦片边长米250 米4800 × 4800231.6563581111950.520500 米2400 × 2400463.3127171111950.5201 公里1200 × 1200926.6254331111950.520于是瓦片原点可以写成X_MIN -20015109.354 Y_MAX 10007554.677 TILE 1111950.520 def tile_origin(h, v): h: 列号 0-35v: 行号 0-17 x_start X_MIN h * TILE y_start Y_MAX - v * TILE return round(x_start, 3), round(y_start, 3) print(tile_origin(0, 0)) # (-20015109.354, 10007554.677) print(tile_origin(35, 17)) # (18903158.846, -8895604.163)拿一块具体瓦片举例列号 25、行号 5 的 1 公里产品左上角 x −20015109.354 25 × 1111950.520 7777243.646左上角 y 10007554.677 − 5 × 1111950.520 4455532.077。把这两个数填进 Map Info 的第四、五项像素大小填 926.6254331投影类型选正弦球体半径填 6371007.181这块瓦片的投影定义就完成了。验证方法是打开图像用光标工具在图像四个角点附近取坐标看数值是否落在上面算出的范围内右下角应该接近 x_start 1200 × 926.625 x_start 1111950.5 左右右上角同理。如果范围对得上投影基本就是对的。如果你习惯在命令行下做交叉验证可以写一小段代码反算经纬度看看落点是不是在那块 10°×10° 的格子内import math R 6371007.181 lon0 0.0 def sinu_to_lonlat(x, y): lat y / R lon x / (R * math.cos(lat)) lon0 return math.degrees(lon), math.degrees(lat) print(sinu_to_lonlat(7777243.646, 4455532.077))算出来的经纬度应该在东经 70° 附近、北纬 40° 附近和瓦片行列号对应的地理位置吻合就说明参数没填错。4. 重投影与 GLT 校正的取舍4.1 GLT 校正适用什么场景正弦投影定义好之后下一步往往是决定要不要转到别的坐标系。这时候有两类做法一类是直接用重投影工具把正弦投影的栅格重采样到地理坐标或投影坐标另一类是先用地理定位文件构建查找表GLT再按查找表重采样。GLT 的思路是先用同一时刻、同一轨道的地理定位产品比如 1 公里的定位数据拿到每个像素对应的经纬度构建一张记录“目标像素该从哪里取值”的索引表然后按这张表去做重采样。它的优势在于利用了原始观测的几何关系对 1B 级这类带几何畸变的数据特别有效边缘和扫描角较大区域的配准精度明显更好。如果是已经做过几何校正的格网产品比如 L3 的反射率、植被指数本身已经在规则正弦格网上不存在几何畸变直接重投影就够了用 GLT 反而是多此一举。判断标准很简单数据是不是规则格网。是的话直接重投影不是且手上有配套定位文件才考虑 GLT。重投影时的重采样方法也要区分对待。分类产品土地覆盖、质量标记必须用最近邻否则会凭空生成不存在的类别值连续型产品反射率、地表温度、植被指数可以用双线性或三次卷积视觉上更平滑统计上也更稳定。4.2 转 WGS84 后位置跑偏的排查路径重投影之后发现图像和矢量边界对不上是最常见的投诉。我一般按下面的顺序排查基本能在十分钟内定位问题。先查源投影定义里球体半径是不是 6371007.181。这一项错了会导致系统性偏移偏移量随纬度变化赤道小、高纬大特征很明显。再查左上角坐标有没有差半像素。ENVI 用的是角点坐标如果你从别的软件里抄来的是像素中心坐标就会整体偏半个像素。然后查重投影时输出范围设置。默认情况下工具会按源数据范围外接矩形生成输出边缘会出现空白带如果后续做了裁剪或掩膜容易误判为位置错误。正确的做法是先看整体范围再看局部。最后查是否混用了不同分辨率的产品。250 米和 1 公里的瓦片虽然覆盖同样的地理范围但像素数差四倍如果定义时把像素大小写反了图像会被缩放位置关系全乱。还有一种隐蔽情况数据本身是经过拼接的全球或大区域图但投影定义还是按单块瓦片的左上角坐标填的。这种错误在小范围看不出来整幅图检查时会发现右侧完全错位。处理办法是重新用全球左上角坐标定义而不是单块瓦片坐标。4.3 常见问题速查表现象最可能原因处理办法图像显示正常但无坐标元数据未解析手工补 Map Info整幅图被横向压扁误定义为经纬度坐标改回正弦投影与矢量叠加整体偏移球体半径填成椭球值改为 6371007.181相邻瓦片拼接有细缝像素大小取整用 231.656358 等精确值重投影后边缘大片空白输出范围按外接矩形按实际有效范围裁剪保存后投影信息消失存成了非 ENVI 格式存为 ENVI 标准格式图层叠加后投影丢失叠加不继承源投影指定参考文件或重定义高纬度区域形状夸张正弦投影固有特性属正常非错误这张表我在实际项目里反复用大部分情况看一眼现象就能对上号。需要强调的是最后一条高纬度形状失真不是错误是正弦投影的数学性质不要试图去“修正”它该做的是换投影。5. 实操心得与几个容易翻车的细节5.1 基准面、球体半径与“差之毫厘”基准面这一项在 ENVI 里容易被当成次要参数忽略。正弦投影用的是自定义球体不是任何一个大地基准面所以在基准面下拉框里应该选“无”或者对应的球体选项而不是顺手选 WGS84。选了 WGS84 的后果是软件在做坐标转换时会按椭球参数去解释那组正弦坐标误差看似不大但在做多源数据融合时会被放大。我做过一次对比同一块瓦片用正确球体半径和用椭球半径分别定义再各自转到地理坐标两者在赤道附近的差异是百米量级到了中纬度扩大到两百米以上。这个量级对于做土地覆盖变化检测、植被指数时序分析来说足以让同一位置的像素落到完全不同的地物上统计结果直接失真。另一个细节是单位。ENVI 的 Map Info 里坐标默认按米解释如果你的数据来源把坐标存成了千米或者其他单位填进去就会差三个数量级。判断方法还是看数值量级正弦投影的全球坐标绝对值在一千万到两千万之间如果看到的是两万左右那基本就是千米单位需要乘 1000。5.2 批处理与结果存档的几点习惯单块瓦片手工定义投影还能接受一批几十上百块的时候就必须脚本化。ENVI 提供了基于 IDL 的批处理能力也可以用它自带的可视化建模工具搭流程。我的做法是把“读元数据、算瓦片原点、写 Map Info、保存”这四步串成一个循环输入是文件列表和瓦片行列号输出是定义好投影的 ENVI 格式文件。这里有个经验批处理时不要直接覆盖原始文件先输出到新目录跑完抽检几块再决定是否清理。投影定义这活儿一旦批量写错回滚成本很高。存档格式也有讲究。定义好投影之后尽量保存为 ENVI 标准格式或者带完整地理信息的 GeoTIFF不要存成普通图像格式。普通格式不保存投影数组等于白干。如果确实需要交付图片单独导出一份原始数据保留完整版本。还有一点关于文件命名。批处理生成的文件建议把瓦片行列号和分辨率编进文件名比如带上 h25v05 和 1km 这样的标识。后期拼全球图或者按区域检索时靠文件名就能定位比打开每个文件看头文件快得多。这个习惯看起来琐碎但数据量上来以后省下的时间非常可观。最后分享一个我常用的自检动作每批数据处理完随机抽一块把图像和一个已知位置的参照数据叠一下看海岸线或者大型水体边界是否吻合。这个动作只要三十秒但能挡住绝大多数投影参数错误。投影这玩意儿定义对了没人夸你定义错了后面所有分析全废花这点时间做验证太值了。