第11章 - 坐标参考系统与高程基准

两批点云叠在一起对不上，X 差了几个数量级；或者叠加得很整齐，但 DEM 算出来的高程整体偏 30 米。前者是水平坐标参考系统（CRS）不匹配，后者是垂直基准不匹配。点云没有底图可以目视校验，一个 SRS 标错要等到土方量、坡度、断面全部算完才会暴露，代价比矢量数据高得多。

本章讲三件事：PDAL 内部如何表示和传递 CRS（SpatialReference 与 PointTable 的规则）；LAS/LAZ 文件里的 SRS 到底从哪些 VLR 读出来、按什么优先级；重投影与高程基准转换时，哪些默认行为会让结果”看起来成功、实际错了”。文中命令与报错均来自 PDAL 2.10.0（本机 git-versioneb9126）实测。

11.1 点云里 SRS 的两个分量

11.1.1 水平 CRS 只是入口

水平 CRS 决定 X/Y 的地球解释。常见两种失误：把投影坐标（单位米、坐标量级 10⁵–10⁷）当成经纬度写进 EPSG:4326；或者把 UTM 33N 的点标成 UTM 32N，横向错开约 500 km 且在局部范围内变形平滑，很难看出来。

PDAL 的做法是把 SRS 作为 Stage 的输出状态之一，而不是附在点坐标上的元数据。一个 stage 处理完 PointView，该 view 会继承 stage 的 SRS（如果有）。这意味着 pipeline 中间某一段把 SRS 改错，后面的 stage 会一路沿用错误 SRS，直到写出文件。

11.1.2 垂直分量才是大头

点云 Z 通常被直接当”高程”用来做 DEM、等高线、土方。但 Z 至少有三种含义：椭球高 h（相对参考椭球，WGS84 下与正高差 −100…+100 m 量级）、正高 H（相对大地水准面）、局部假定高程（施工坐标，可能任意偏移）。把椭球高当正高用，在中国东部、欧洲、北美都有几十米级系统偏差。

PDAL 没有单独的”垂直基准”stage。垂直信息只能靠两处表达：SRS 本身是复合 CRS（例如 EPSG:4326+5773 表示 WGS84 水平 + EGM96 高程），或者交给 PROJ 的 pipeline/网格参数。第 11.5 节会用一个实测例子说明：复合 CRS 写对了，垂直分量也可能被静默忽略。

11.1.3 SRS 在 pipeline 中的流转

reader (LAS VLR / SRS 选项)
   │  SpatialReference
   ▼
Stage::setSpatialReference ──► PointView 继承 stage 的 SRS
   │
   ├─ filters.reprojection  改 X/Y/Z，并把 stage SRS 改成 out_srs
   ├─ filters.ferry         只搬维度，不动 SRS
   ├─ filters.assign        只改属性，不动 SRS
   ▼
writer (a_srs 优先，否则用输入 view 的 SRS 写 VLR)

图中的规则不是约定，而是实现：LasWriter::prepOutput() 里就是 m_srs = getSpatialReference().empty() ? srs : getSpatialReference()——writer 自己设了 a_srs 就用它，否则用上游 view 的 SRS。

11.2 SpatialReference 的解析规则

11.2.1 字符串怎么变成 SRS

pdal::SpatialReference::set() 对输入分三路处理：

第三种是绝大多数情况。解析失败抛出的错误长这样：

PDAL: Could not import coordinate system 'EPSG:99999': PROJ: proj_create_from_database: crs not found

注意它来自 GDAL 的 SetFromUserInput，所以 GDAL/PROJ 能认的写法都能用；反过来，PROJ 数据库里没有的 EPSG 代码一定失败，和 PDAL 版本无关。

11.2.2 一个点表只能有一个”确定”的 SRS

SpatialReference.hpp 的类注释写得很清楚：点表跟踪当前处理中的 view 的 SRS；如果一个点表里有多个不同的 SRS，PointTable::spatialReference() 返回空 SRS，spatialReferenceUnique() 返回 false。

这直接决定了一个高频故障：用 filters.merge 合并两个 SRS 不同的文件，merge 后点表 SRS 变空，紧接着的 filters.reprojection 要么报错，要么被迫用 in_srs 强行指定（那一半数据必然错）。正确做法是先逐文件重投影到同一 out_srs，再合并。

11.2.3 查看一个文件真实的 SRS

pdal info --summary 会输出 metadata.srs 子树，这是核对 SRS 最省事的方式：

pdal info utm33.las --summary

关键字段（PDAL 2.10.0 实测）：

"metadata": {
  "srs": {
    "horizontal": "PROJCS[\"WGS 84 / UTM zone 33N\", ...]",
    "vertical": "",
    "compoundwkt": "PROJCS[\"WGS 84 / UTM zone 33N\", ...]",
    "proj4": "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs",
    "isgeographic": false,
    "isgeocentric": false,
    "units": { "horizontal": "metre", "vertical": "" },
    "json": { "type": "ProjectedCRS", "id": { "authority": "EPSG", "code": 32633 }, ... }
  },
  "spatialreference": "PROJCS[\"WGS 84 / UTM zone 33N\", ...]",
  "comp_spatialreference": "PROJCS[...]"
}

判读规则：vertical 为空字符串就说明这个文件没有垂直基准声明，Z 的来源是未知的（多半是椭球高或仪器高）。units.vertical 同理。json 是 PROJJSON，比 WKT 好解析，做批量 SRS 审计时优先用它。

11.3 LAS/LAZ 里的 SRS 从哪读

11.3.1 三套 VLR 与优先级

LAS 文件里 SRS 可以存在三套地方：GeoTIFF keys（LASF_Projection 下的 34735/34736/34737 号 VLR）、WKT1（34735 之外的 WKT record，LAS 1.4 强制）、以及 PDAL 专用的 WKT2/PROJJSON VLR。io/private/las/Srs.cpp 的 Srs::init() 规定默认顺序：

mustUseWkt() 的定义在 io/private/las/Header.hpp：WKT global-encoding 位被置且 minor version ≥ 4，或点格式 > 5。也就是说 LAS 1.2 + 点格式 3 的文件走 GeoTIFF 兜底，LAS 1.4 走 WKT 兜底。

先命中的赢，后续不再尝试。这在混用软件产出的文件上会出问题：某些工具同时写了 GeoTIFF 和 WKT，两者不一致时 PDAL 只取一个，且不报错（Debug 日志里只有一句 “File contains both WKT and GeoTiff VLRs which is disallowed.”）。

11.3.2 用 srs_vlr_order 强制选择

readers.las 和 readers.copc 都有 srs_vlr_order 选项，接受逗号分隔的 wkt2,wkt1,projjson,geotiff。非法值直接报错：

PDAL: readers.las: Invalid SRS type 'wkt3'. Must be one of 'wkt1', 'geotiff', 'projjson', 'wkt2'.

如果指定了顺序但一个有效 VLR 都没找到，只有一条 Warning：'srs_vlr_order' specified but no valid VLR was found.，pipeline 继续跑，只是数据没有 SRS。批量处理时这类 warning 很容易被日志噪音淹没，建议把 --loglevel 调到 info 并抓关键词。

11.3.3 写出时如何决定 VLR

writers.las 的行为可以归结成三条：

1. a_srs 非空 → 用 a_srs（initialize() 里 setSpatialReference(aSrs)）；
2. 否则用上游 view 的 SRS；
3. LAS 1.4 输出总是置 WKT global-encoding 位并写 WKT1 记录；enhanced_srs_vlrs 打开（默认 false）时额外写 PDAL 私有 WKT2 与 PROJJSON VLR。

还有一个容易误解的点：forward 不转发 LASF_Projection 容器 VLR。源码注释写得很明确——这些 VLR 由 writer 自己重建，所以你不能靠 forward="vlr" 把源文件的 GeoTIFF keys 原样搬到输出。

最小可跑示例，把一个无 SRS 的合成文件标定并转成 LAS 1.4：

[
  {
    "type": "readers.faux",
    "filename": "faux",
    "count": 2000,
    "mode": "Normal",
    "mean_x": 400250, "mean_y": 5700250, "mean_z": 35,
    "stdev_x": 80, "stdev_y": 80, "stdev_z": 10,
    "override_srs": "EPSG:32633"
  },
  {
    "type": "writers.las",
    "filename": "utm33.las",
    "minor_version": 4,
    "dataformat_id": 7,
    "scale_x": 0.01, "scale_y": 0.01, "scale_z": 0.01,
    "offset_x": "auto", "offset_y": "auto", "offset_z": "auto",
    "enhanced_srs_vlrs": true
  }
]

readers.faux 不是玩具才有用：写回归测试时它是唯一不依赖外部数据的选择。上面 mode 只接受 Ramp、Uniform、Normal、Grid 四种；写错会得到 Option 'seed' not supported with mode 'RampUniformNormalGridInvalid'.——报错里那串枚举就是源码里的全部合法值。

11.4 filters.reprojection 的真实行为

11.4.1 机制与必填项

filters.reprojection 用 GDAL 的 OGRCoordinateTransformation 逐点变换 X/Y/Z，旧坐标被直接替换。要保留原坐标，先用 filters.ferry 复制维度：

[
  "utm33.las",
  { "type": "filters.ferry", "dimensions": "X=>Xutm, Y=>Yutm" },
  { "type": "filters.reprojection", "out_srs": "EPSG:4326" },
  { "type": "writers.las", "filename": "wgs84.las",
    "scale_x": 0.0000001, "scale_y": 0.0000001, "scale_z": 0.01,
    "offset_x": "auto", "offset_y": "auto", "offset_z": "auto" }
]

参数（源码 filters/ReprojectionFilter.cpp 的 addArgs()）：

11.4.2 输入没 SRS 时的报错

如果输入 view 没有 SRS，又没有给 in_srs，报错文本是固定的（实测）：

PDAL: filters.reprojection: source data has no spatial reference and none is specified with the 'in_srs' option.

修法按优先级：最好的做法是源头修正（reader 用 override_srs 或 default_srs 声明真实 SRS）；其次是 pipeline 里显式给 in_srs；两者都别做”猜一个”——猜错的投影不会有任何提示。

override_srs 与 default_srs 是所有 reader 的通用选项（doc/stages/reader_opts.md），二者不能同时使用：前者无条件覆盖文件里的 SRS，后者只在文件没有 SRS 时补上。

11.4.3 坐标精度：默认 scale 会把地理坐标压扁

writers.las 的默认 scale_x/y/z 都是 0.01。这在投影坐标（米）下是 1 cm，在地理坐标（度）下是 0.01° ≈ 1.1 km。实测把一个 2000 点的 UTM 33N 文件转到 EPSG:4326 且不指定 scale：

pdal info wgs84_defaultscale.las --summary
# "minx": 13.56, "maxx": 13.57,
# "miny": 51.44, "maxy": 51.45,
# "scale_x": 0.01, "scale_y": 0.01,

原始范围在 500 m × 500 m 量级，转成度后只剩 0.01° 的格网——所有点被量化到几十个格子上，几何信息基本丢光，而 pipeline 不报任何错误。转到地理坐标时 scale_x/scale_y 应给 0.0000001（约 1 cm）。

11.4.4 轴序：PDAL 默认与传统 GIS 一致

GDAL 3 起，WKT 里的轴序会被严格执行，导致 EPSG:4326 按 lat/lon 而非 lon/lat 解释。PDAL 在 SrsTransform::set() 里对源和目标都调用 SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER)，所以默认就是 X=经度、Y=纬度，不需要像调 GDAL 命令行那样额外加 -s_srs 轴序翻转。

只有在数据本身存的就是 lat/lon 时，才需要 in_axis_ordering：

{ "type": "filters.reprojection", "in_srs": "EPSG:4326", "out_srs": "EPSG:32633",
  "in_axis_ordering": ["2", "1"] }

轴序向量的三个校验都来自源码：元素必须能转成整数（否则 Unable to convert axis ordering to integer）、长度 ≤ 3、排序后必须是 1,2 或 1,2,3。

11.4.5 失败点会被静默丢弃

error_on_failure 默认 false。逐点变换失败（超出投影有效范围、反投影无解）时，processOne() 返回 false，该点不进入输出 view，pipeline 成功退出。表现是输出点数比输入少，且 summary 里只有 count 一个数字能看出来。

生产 pipeline 里我建议一律打开：

{ "type": "filters.reprojection", "out_srs": "EPSG:4526", "error_on_failure": true }

打开后的错误带坐标，便于定位脏数据：

PDAL: filters.reprojection: Couldn't reproject point with X/Y/Z coordinates of (nan, nan, nan).

11.5 高程基准：PDAL 不做转换，PROJ 做

11.5.1 没有 vertical datum stage

PDAL 的阶段里没有 filters.verticaldatum 之类的东西。垂直基准转换只有两条路：

1. 复合 CRS：out_srs 写成 EPSG:4326+5773（EGM96）或国家复合代码（如丹麦 EPSG:7416、中国 2000 国家高程基准对应的地方代码），由 PROJ 查库决定变换；
2. 显式 pipeline：filters.projpipeline 的 coord_op 直接给 PROJ 串，例如 +proj=pipeline ... +proj=vgridshift +grids=egm96_15.tif，或 +proj=affine +zoff=100。

filters.projpipeline 的三个参数：

注意 out_srs 那一行：文档明确说”No check is done to ensure the compliance”。你可以把一个米制平移写成 +proj=affine +zoff=-32.5，同时声明 out_srs 是某个正高复合 CRS，PDAL 照写不误。这既是灵活性，也是最容易产出”元数据撒谎”文件的入口。

该 stage 依赖 GDAL ≥ 3.0，2.x 的 GDAL 环境下二进制里根本没有这个 driver。

11.5.2 一个代价很高的实测

环境：PROJ 9.8.1（OSGeo4W），PROJ_NETWORK=OFF，本机没有 EGM96 网格。把 UTM 33N 的椭球高数据转成 EPSG:4326+5773：

import pyproj
from pyproj import CRS, Transformer
import pyproj.network
pyproj.network.set_network_enabled(False)

src = CRS.from_epsg(32633)
dst = CRS.from_user_input("EPSG:4326+5773")   # WGS84 + EGM96 height
t = Transformer.from_crs(src, dst, always_xy=True)
print(t.description)
print(t.transform(400000.0, 5700000.0, 50.0))

输出：

Inverse of UTM zone 33N + Inverse of Null geographic offset from WGS 84 to WGS 84 (with axis order normalized for visualization)
(13.561123205604062, 51.44235231637489, 50.0)

变换描述里的 Null geographic offset 是关键：PROJ 找到了复合 CRS 的定义，但缺少 EGM96 网格文件，于是垂直分量退化为”空偏移”，Z 原封不动返回 50.0。整个过程没有异常、没有 warning，pipeline 会正常写出一个自称是 EGM96 正高、实际是椭球高的文件。

判断方法只有两个：看 projinfo 给出的 pipeline 里有没有 vgridshift 步骤，或看 pdal info --summary 的 metadata.srs.vertical 是否非空（非空只证明声明了垂直分量，不证明网格存在）。

# 看 PROJ 实际会走哪条 pipeline
projinfo -s EPSG:32633 -t EPSG:4326+5773 --spatial-test intersects

# 看 PROJ 的数据目录和网格搜索路径
projinfo --searchpaths

11.5.3 网格文件的部署

PROJ 9 之后网格是 GeoTIFF（旧格式 .gtx 仍被支持）。三种部署方式：

PDAL 仓库自带的 doc/tutorial/nz_reproject.json 与 grid-shift.md 用的就是第三种：把新西兰本地高程基准的 GTX 文件写进 filters.reprojection 的 out_srs 字符串，再在 writers.las 里用 a_srs 声明最终复合代码，否则 PDAL 会把垂直基准写成 “Unnamed Vertical Datum”（文档原话）。这一条值得照抄——先转换、再声明，两步分开做。

11.6 排错清单

11.7 生产建议

先审计再处理。 一批数据入库前跑一遍 pdal info --summary，把 metadata.srs.json.id（PROJJSON 里的 authority/code）和 units.vertical 抽出来列表。同一次任务里出现两种水平代码就直接停下，不要用 in_srs 硬转——那等于把错配写进结果。

中间结果保持投影坐标系。 米制坐标下 1 cm 精度只需 scale 0.01，地理坐标下需要 1e-7；而且面积、坡度、距离算法在投影坐标上才成立。只有最终交付或入库 PostGIS geography 时才转经纬度。

垂直基准单独作为一条处理链。 不要把”水平重投影”和”椭球高→正高”塞进同一个 stage 配置里；先确认网格文件到位（projinfo 看得到 vgridshift），再做垂直转换，最后用 a_srs 声明结果。三者分开，出错时能定位到具体一步。

给 reprojection 打开 error_on_failure。 这是唯一能让”丢点”变成”报错”的开关，成本是失败即中断——批处理脚本把它包在 try/except 里逐个文件重试即可。

LAS 1.4 优先。 1.4 强制 WKT VLR，能表达复合 CRS 和三维基准；1.2 走 GeoTIFF keys，很多复合信息根本存不下。代价是文件头和每个点的结构都比 1.2 大，但对需要 SRS 正确性的工程可以接受。

11.8 本章小结