PostGIS 里处理栅格的人,迟早会遇到一个特别磨人的问题:两张栅格单独看都正常,波段范围对、NoData 也设置了、值域也合理,但只要放到一起按像元做计算,结果就莫名其妙地偏。我第一次踩这个坑是在做两期影像的变化检测,差值图出来成了一张“条纹码”,后来才发现根本不是算法写错了,而是两张栅格的网格根本没有对齐。PostGIS 专门提供了一个函数叫ST_SameAlignment,就是用来判断两个栅格是否落在同一套网格上的。它本身不复杂,但在我接触到的很多项目里,真正理解它判断逻辑的人并不多。这篇东西就把对齐检查的用途、判断细节、实际写法和踩过的坑一次讲清楚,适合准备做栅格叠加分析或者正在被错位结果折磨的朋友。
1. 为什么要优先处理栅格对齐问题
1.1 栅格对齐在 PostGIS 里的具体含义
栅格本质上就是一张坐标纸,每个像元有自己的横纵坐标范围。所谓对齐,不是看两张图分辨率一样,也不是表面上看范围差不多,而是要求两张坐标纸的每一个格子都能严丝合缝地叠在一起。PostGIS 里的栅格对象,除了保存像素值矩阵之外,还会保存一套空间参考信息,包括左上角原点坐标、横纵像元大小、倾斜参数、SRID 等等。只有当这些元数据所确定的网格完全重合时,两个栅格才能按像元一一对应。
这个网格重合的判断标准,比很多人想象的要严格。两张栅格分辨率完全一样,但左上角水平方向差了半个像元,那它们对齐吗?不对齐。为什么?因为每个格子的边界都错开了,你在 A 的某个像元里取值,对应到 B 已经是两个像元中间的位置。反过来,如果两张栅格的左上角不一样,但偏移量恰好是像元尺寸的整数倍,那它们依然是对齐的。这一点在后面实操部分会再用 SQL 验证。
1.2 不对齐对后续运算的直接影响
不对齐最直接的后果,就是像元级运算结果失真。比如计算差值、计算比值、做分类统计,PostGIS 是在每个像元位置上取两个栅格的值做计算。如果两个栅格的格子没有对齐,那么参与计算的并不是同一地表位置,差值自然没有任何物理意义。
另一个容易被忽略的问题是重采样。部分函数在发现两个栅格不对齐的时候并不会直接报错,而是内部悄悄做一次重采样,把 B 拉到 A 的网格上。听起来好像很智能,但这种隐式重采样的算法、每次输入顺序不同导致的差异,经常让人防不胜防。你以为自己在处理原始数据,实际上在比较一个已经被插值处理过的结果和另一个未被处理的原始结果,隐患非常大。所以我在项目里定了一条规矩:任何双栅格运算之前,必须先用ST_SameAlignment亮明身份,没对齐就显式处理,绝不把重采样决定权交给函数内部的默认行为。
1.3 ST_SameAlignment 的定位
ST_SameAlignment的作用就是快速回答一个布尔问题:这两个栅格是不是同一个网格。它不需要你手动比较七八个元数据字段,也不用你编写带容差的判断逻辑,一个函数调用就搞定。虽然你完全可以用ST_ScaleX、ST_ScaleY、ST_SkewX、ST_SRID这些函数自己拼一个判断,但在代码可读性和维护性上,原生函数显然更合适。
不过它也只是一个“体检指标”,只回答网格是否对齐,不回答数据质量、波段结构、NoData 是否一致。实际使用中,我通常把它放在检查流程的第一步:先看网格,再看波段,再看取值范围。网格没对齐,后面所有比较都免谈。
2. ST_SameAlignment 的核心判据与注意事项
2.1 函数形式与一条最简单的测试 SQL
函数签名很直接:
boolean ST_SameAlignment(raster rastA, raster rastB)返回true表示两个栅格对齐,返回false表示不对齐。下面这条 SQL 是全篇最核心的用法,也最适合拿来感受函数行为:
SELECT ST_SameAlignment( ST_MakeEmptyRaster(10, 10, 500000, 4000000, 30, -30, 0, 0, 32650), ST_MakeEmptyRaster(10, 10, 500000, 4000000, 30, -30, 0, 0, 32650) ) AS aligned;ST_MakeEmptyRaster的参数依次是宽度、高度、左上角 X、左上角 Y、像元宽度、像元高度、X 方向倾斜、Y 方向倾斜、SRID。这时候两张栅格所有参数完全一致,查询结果自然是true。
2.2 它检查的五个要素
ST_SameAlignment的判断可以拆成这么几块:
| 检查要素 | 含义 | 不满足时的结果 |
|---|---|---|
| 像元宽度 | 每个格子在 X 方向的实际尺寸 | 一个栅格被横向拉伸或压缩,不能对应 |
| 像元高度 | 每个格子在 Y 方向的实际尺寸 | 纵向对应关系错位 |
| 倾斜参数 | 栅格是否发生旋转、剪切 | 旋转后的栅格与原栅格网格交叉 |
| SRID | 空间参考系 | 坐标系不同,坐标值无法直接比较 |
| 网格原点 | 两个左上角原点是否在同一个格网上 | 即使分辨率相同也可能相位错位 |
其中最容易忽略的是“网格原点”这一项。ST_SameAlignment并不要求两张栅格的左上角完全相等,而是要求两个左上角坐标的差值,能够被像元尺寸整除。换句话说,只要一个栅格相对另一个在 X、Y 方向平移了整数个像元,它们依然能落在同一个无限格网上,这种数据也是对齐的。
对于 Y 轴还要特别留意,PostGIS 栅格的像元高度经常是负数,表示从左上角开始向下递增,所以在手写公式时最好取绝对值。例如:
dx / |scaleX| 必须是整数 dy / |scaleY| 必须是整数skew 参与的情况下判断会更复杂,但实际生产数据多数都是 0。如果你处理的数据带地理参考旋转,直接信任ST_SameAlignment就好。
2.3 它不检查什么
ST_SameAlignment只关心几何网格,不关心你存了什么内容。下面这些因素不会影响它的判断结果:
- 波段数量和波段顺序
- 像素类型,比如
8BUI还是32BF - NoData 值设置
- 像素值的实际大小
- 栅格的宽度和高度
- 两个栅格是否真的有空间重叠
最后一条尤其容易造成误解。A 和 B 分别位于两个互不相邻的区域,但它们的网格定义完全一致,ST_SameAlignment依然返回true。对齐是对几何坐标系的描述,不是对相交区域的描述。这也提醒我们:对齐检查通过了,下一步仍然要判断空间范围关系是不是符合你的运算目标。
2.4 边界情况:网格重合但位置不同
我举个例子。三个栅格都使用 30 米分辨率,X 方向有效范围大概是 500000 到 500300 左右,但左上角 X 坐标分别是 500000、500060、500075:
WITH r AS ( SELECT 1 AS rid, ST_MakeEmptyRaster(10, 10, 500000, 4000000, 30, -30, 0, 0, 32650) AS rast UNION ALL SELECT 2 AS rid, ST_MakeEmptyRaster(10, 10, 500060, 3999940, 30, -30, 0, 0, 32650) AS rast UNION ALL SELECT 3 AS rid, ST_MakeEmptyRaster(10, 10, 500075, 3999940, 30, -30, 0, 0, 32650) AS rast ) SELECT a.rid, b.rid, ST_SameAlignment(a.rast, b.rast) AS aligned FROM r a JOIN r b ON a.rid < b.rid ORDER BY a.rid, b.rid;结果会是这样:
| a.rid | b.rid | aligned |
|---|---|---|
| 1 | 2 | true |
| 1 | 3 | false |
| 2 | 3 | false |
1 和 2 的 X 差 60 米,正好是 2 个像元;Y 差 60 米,也正好是 2 个像元,所以对齐。1 和 3 的 X 差 75 米,等于 2.5 个像元,网格相位不对,返回false。这组数据可以留着以后做实验,帮助自己直观理解“相位错位”是怎么回事。
3. 实操:用 ST_SameAlignment 做数据体检
3.1 准备对齐测试数据
如果你的数据库里还没有现成栅格表,最简单的办法是用ST_MakeEmptyRaster创建元数据栅格,再用ST_AddBand加一个波段。比如:
CREATE TABLE test_rasters ( rid integer PRIMARY KEY, rast raster ); INSERT INTO test_rasters VALUES (1, ST_AddBand(ST_MakeEmptyRaster(10, 10, 500000, 4000000, 30, -30, 0, 0, 32650), '8BUI'::text, 0, 0)), (2, ST_AddBand(ST_MakeEmptyRaster(10, 10, 500060, 3999940, 30, -30, 0, 0, 32650), '8BUI'::text, 0, 0)), (3, ST_AddBand(ST_MakeEmptyRaster(10, 10, 500075, 3999940, 30, -30, 0, 0, 32650), '8BUI'::text, 0, 0));实际生产数据一般来自raster2pgsql,但元数据结构完全一致。测试表能帮你快速理解函数行为,方便验证后面所有 SQL。
3.2 查单对栅格的对齐状态
数据就绪之后,最简单的检查:
SELECT ST_SameAlignment( (SELECT rast FROM test_rasters WHERE rid = 1), (SELECT rast FROM test_rasters WHERE rid = 2) ) AS aligned_result;返回true。但仅知道结果还不够,调试时我更推荐把关键的元数据字段全部打出来。这样一旦结果不符合预期,可以立刻看出是哪一项出了问题:
SELECT ST_SRID(a.rast) AS srid_a, ST_SRID(b.rast) AS srid_b, ST_ScaleX(a.rast) AS scalex_a, ST_ScaleX(b.rast) AS scalex_b, ST_ScaleY(a.rast) AS scaley_a, ST_ScaleY(b.rast) AS scaley_b, ST_SkewX(a.rast) AS skewx_a, ST_SkewX(b.rast) AS skewx_b, ST_UpperLeftX(a.rast) AS ulx_a, ST_UpperLeftX(b.rast) AS ulx_b, ST_UpperLeftY(a.rast) AS uly_a, ST_UpperLeftY(b.rast) AS uly_b, ST_SameAlignment(a.rast, b.rast) AS is_aligned FROM test_rasters a, test_rasters b WHERE a.rid = 1 AND b.rid = 2;这种“ST_SameAlignment 给结论,手工字段给证据”的组合,是我排查栅格问题最常用的方式。
3.3 批量检查一张瓦片表
单对检查满足不了真实项目。实际生产环境里可能有一整张瓦片表,几十万块瓦片,你很难手动挑几对检查。
当表规模不大,比如几百块瓦片时,可以暴力检查所有需要参与运算的相邻瓦片对:
SELECT a.rid AS rid_a, b.rid AS rid_b, ST_SameAlignment(a.rast, b.rast) AS aligned FROM test_rasters a JOIN test_rasters b ON b.rid > a.rid AND ST_Intersects(ST_Envelope(a.rast), ST_Envelope(b.rast)) ORDER BY a.rid, b.rid;ST_Envelope把栅格的外包矩形转成几何对象,再用ST_Intersects限制到相邻瓦片,避免出现完全不相关的瓦片两两比较。这个查询在几十块的测试表上完全够用,但如果瓦片数量上百,建议换一种思路:不查全表,只抽样本或者只查边界互相接触的瓦片。
批量检查结果我最想看到的是四种情况:全对齐,少数不对齐,成片不对齐,以及混合 SRID。每一种的后续处理策略不一样。全对齐直接进入计算;少数不对齐优先检查是不是浮点误差;成片不对齐通常是数据源头或投影方式不统一;混合 SRID 则必须先做坐标系统一。
3.4 把检查结果放进地图代数流程
ST_SameAlignment最实用的场景,是作为地图代数运算的前置条件。下面是一个简单例子,对两张对齐栅格做像元均值:
UPDATE target_result AS t SET rast = ST_MapAlgebraExpr(a.rast, 1, b.rast, 1, '([rast1] + [rast2]) / 2.0', '32BF'::text) FROM test_rasters a, test_rasters b WHERE a.rid = 1 AND b.rid = 2 AND ST_SameAlignment(a.rast, b.rast);ST_MapAlgebraExpr的表达式里,[rast1]对应第一个栅格,[rast2]对应第二个栅格。如果ST_SameAlignment返回false,这条UPDATE就不会执行,至少不会在你没准备的情况下发生隐式重采样。这种防御式写法虽然不能解决不对齐问题,但能让失败提前暴露,不会让脏结果溜进下一层。
3.5 不对齐时的补救:重采样到参考栅格
一旦发现不对齐,我一般的策略是选一张“基准栅格”,然后把所有其他栅格通过ST_Resample转换到基准网格上。
UPDATE target_tiles AS t SET rast = ST_Resample( ST_Transform(src.rast, ST_SRID(t.rast)), t.rast, 'NearestNeighbor' ) FROM source_tiles AS src WHERE src.rid = 1 AND t.rid = 100;先ST_Transform把坐标系统一到和目标栅格一致,再用ST_Resample以目标栅格t.rast作为参考网格,把源栅格重新采样进去。第三个参数NearestNeighbor是重采样算法。分类数据、土地覆盖数据一定要用最近邻,因为双线性或三次卷积会插出根本不存在的类别值;连续变量比如温度、高程、NDVI,可以用Bilinear。执行完这张更新之后,不要直接相信结果,再跑一次ST_SameAlignment确认一下,这是我最常提醒自己的事。
4. 常见坑与排查技巧
4.1 坐标系转换后“必然不对齐”的坑
最典型的问题是ST_Transform之后直接参与计算。把栅格从 UTM 投影转成 Web Mercator,或者从地理坐标转成平面坐标,每个像元在目标坐标系下的宽度已经不是整齐的整数米了。即便你在 SQL 里指定了输出分辨率,重投影本身很可能产生亚像素偏移。
我曾经处理过一次数据:源数据是 UTM 50N,10 米分辨率,目标表是 UTM 50N,也是 10 米分辨率。看起来同一个投影,应该没问题。但源数据经过了ST_Transform,输出的像元原点比目标表多了 0.000000001 度级别的残差。单看没有任何问题,一算ST_SameAlignment就是false。后来我的标准流程改成:先ST_Transform,再ST_Resample到目标网格,最后再用ST_SameAlignment校验,三步缺一不可。
有一个细节值得记住:ST_Transform只负责坐标换算,不负责让目标栅格和现有栅格对齐。很多网上教程只提到ST_Transform用来转坐标系,却忘了重投影之后必须重采样,这就是很多人踩坑的原因。
4.2 浮点精度造成的假失败
浮点误差是ST_SameAlignment返回false的另一个高频原因。比如两张栅格名义上都是 30 米分辨率,但一个存的是29.999999999,另一个存的是30.000000001,在 PostGIS 内部用双精度比较时可能就会判为不对齐。
遇到这种情况,不要直接放弃。先跑参数诊断 SQL,看每一项差了多少。如果只是 1e-9 级别的差异,我建议写一个带容差的辅助函数,而不是死磕ST_SameAlignment的严格返回值。下面这个函数是我在项目里常用的版本:
CREATE OR REPLACE FUNCTION raster_grid_aligned( r1 raster, r2 raster, eps double precision DEFAULT 1e-6 ) RETURNS boolean LANGUAGE sql IMMUTABLE AS $$ SELECT ST_SRID(r1) = ST_SRID(r2) AND abs(ST_ScaleX(r1) - ST_ScaleX(r2)) < eps AND abs(ST_ScaleY(r1) - ST_ScaleY(r2)) < eps AND abs(ST_SkewX(r1) - ST_SkewX(r2)) < eps AND abs(ST_SkewY(r1) - ST_SkewY(r2)) < eps AND abs( (ST_UpperLeftX(r1) - ST_UpperLeftX(r2)) / abs(ST_ScaleX(r1)) - round((ST_UpperLeftX(r1) - ST_UpperLeftX(r2)) / abs(ST_ScaleX(r1))) ) < eps AND abs( (ST_UpperLeftY(r1) - ST_UpperLeftY(r2)) / abs(ST_ScaleY(r1)) - round((ST_UpperLeftY(r1) - ST_UpperLeftY(r2)) / abs(ST_ScaleY(r1))) ) < eps; $$;eps默认取 1e-6,表示允许 0.000001 个像元的偏差。这个容差不要放大太多,否则真的不对齐的数据也会被放过去。这里的容差单位是像元个数,不是实际距离,所以跟分辨率无关,这也是我推荐这种写法的原因。
4.3 “分辨率相同但相位完全不同”的坑
还有一类数据更隐蔽:分辨率一样,SRID 一样,甚至左上角坐标看起来很接近,但差了几十米,落不到同一格网上。比如两张遥感影像,一个左上角是 (500020, 4000000),另一个是 (500045, 3999985),都是 30 米分辨率。前者偏移了 20 米,后者偏移了 45 米,两个都不是 30 的整数倍,网格完全错开。这种情况不会报“投影错误”,也不会影响单张图显示,但一旦叠加计算,结果就会花掉。
我处理这种数据时有一个习惯:先看差异距离是否小于一个像元。如果差异本身小于像元尺寸,说明是数据生产过程中的裁剪偏移,可以考虑用ST_Resample统一到某个网格;如果差异是几十个像元,那你得搞清楚是不是某张图用了不同的切片方案。网格不是你想当然的那样,切瓦片的原点设置、行列数、以及是否从某个固定起点开始,都会影响原点位置。做全球或者全国遥感数据时,先查一下生产规范里的格网原点定义,能少走很多弯路。
还有一个额外建议:在做像元对齐之前,先把ST_ScaleX和ST_ScaleY检查完。只要分辨率不一致,你再怎么平移原点也没用,只能重采样。所以排查顺序永远是“SRID → 分辨率 → 倾斜 → 原点”。
4.4 不同 SRID 的同网格数据
有两张影像,一张用 EPSG:32650,另一张用 EPSG:32651。虽然两个坐标系的中央经线不同,但在相邻区域它们的数值会差很大,直接比较没有任何意义。ST_SameAlignment看到 SRID 不同,直接返回false。
这类问题其实最好办,做一次ST_Transform统一到同一个 SRID 就行。但要注意,转换后的栅格很可能直接变成不对齐状态,所以仍然要走一遍ST_Resample。不要以为 SRID 统一了就万事大吉,统一 SRID 只是第一步,网格对齐才是第二步。我在实际项目中见过不少人只做了ST_Transform就跑地图代数,最后结果照样不对。
4.5 空瓦片与 NoData 对检查结果的影响
栅格瓦片有可能会遇到完全没有任何有效像元的空瓦片。空瓦片参与ST_SameAlignment可能产生三种状态:true、false,甚至NULL。如果返回NULL,直接进WHERE判断时会被当成假值,容易让报告里出现“全表没有几对对齐”的假象。
所以跑批量检查之前,建议先过滤掉空瓦片。PostGIS 里可以用ST_Count(rast)统计有效像元数,只保留大于 0 的瓦片:
WHERE ST_Count(a.rast) > 0 AND ST_Count(b.rast) > 0NoData 对ST_SameAlignment本身没有影响,但对后续计算影响非常大。两张栅格对齐了,但一张的 NoData 区域和另一张的有效区域重叠,运算结果里会出现大片异常值。我的经验是:先对齐,再统一 NoData 值,最后再计算。顺序反了,排查成本会成倍增加。
5. 我的几点实操体会
5.1 先在入库前统一网格,而不是事后补救
PostGIS 栅格处理有点像盖房子,地基没打好,后面所有装修都白费。我目前比较推荐的方式是在数据写进数据库之前,先用 GDAL、QGIS 或者其他工具把分幅影像统一重采样到同一个网格规范。入库之后,在每个瓦片插入时顺手检查一遍ST_SameAlignment,不对齐的直接拦下来,不要等运算阶段再发现。
这样做看起来多了一步,实际上能省掉后面大量调试时间。有一次我在处理一个区域的 DEM 拼接,由于各源文件分辨率有 1 米和 2 米两种,我直接在入库脚本里加入对齐校验,凡是没对齐的自动转成 2 米网格。整个流程跑完只花了几分钟,但避免了后期写一堆重采样逻辑的麻烦。
5.2 检查报告要保留到数据交付
在正式项目里,我不仅会跑ST_SameAlignment检查,还会把检查结果落成一张表。哪两张瓦片对齐、哪些瓦片经过了重采样、用了什么算法、容差是多少,全部记录在案。数据交付时,这张表本身就是质量报告的一部分。别人拿到你的数据可以先看报告,而不是重新发现问题再来找你。
这类报告用 SQL 写起来其实很简单,核心思路是把参数差异和检查结果输出成行:
SELECT a.rid AS rid_a, b.rid AS rid_b, ST_SRID(a.rast) AS srid_a, ST_SRID(b.rast) AS srid_b, ST_ScaleX(a.rast) AS scalex_a, ST_ScaleX(b.rast) AS scalex_b, ST_UpperLeftX(a.rast) - ST_UpperLeftX(b.rast) AS ulx_diff, ST_UpperLeftY(a.rast) - ST_UpperLeftY(b.rast) AS uly_diff, ST_SameAlignment(a.rast, b.rast) AS aligned FROM tiles a JOIN tiles b ON b.rid > a.rid WHERE ST_Count(a.rast) > 0 AND ST_Count(b.rast) > 0 ORDER BY aligned, a.rid, b.rid;这里的ulx_diff和uly_diff可以直接看出原点差是否接近像元尺寸的整数倍。
5.3 最后留一个快速巡检 SQL
最后分享一个我在每次数据处理完必跑的快速巡检 SQL。它不替代ST_SameAlignment,但能从全局角度快速发现异常瓦片:
SELECT ST_SRID(rast) AS srid, round(ST_ScaleX(rast)::numeric, 6) AS scalex, round(ST_ScaleY(rast)::numeric, 6) AS scaley, count(*) AS tile_count, count(DISTINCT round(ST_UpperLeftX(rast)::numeric, 6)) AS distinct_ulx, count(DISTINCT round(ST_UpperLeftY(rast)::numeric, 6)) AS distinct_uly FROM tiles GROUP BY 1, 2, 3 ORDER BY 1, 2, 3;如果distinct_ulx或distinct_uly数量很多,说明原点五花八门,十有八九存在对齐问题。这种聚合查询非常轻量,适合放进定时任务,每次新增数据后跑一遍。栅格对齐这件事不复杂,但它是一个典型的“小问题酿大事故”的环节。先把ST_SameAlignment用熟、用对,再配合容差判断和重采样流程,你处理 PostGIS 栅格数据的效率会明显提升,至少不会在最基础的网格上浪费时间。