第08章 - 几何裁剪与重投影过滤器

8.1 两类最常被混用的 filter

裁剪（crop）与重投影（reprojection）在流程里经常挨着出现，但它们的性质完全不同：

• filters.crop 是丢点操作。它不改变任何坐标，只是决定哪些点留下；
• filters.reprojection 是改坐标操作。它一个点都不丢，但会重写 X/Y/Z。

所以它们在 pipeline 里的位置很明确：先裁剪后重投影。裁剪只是对原坐标做比较，重投影要对每个点做一次变换（含 datum 变换时更贵）；而且裁剪范围通常在原始坐标系下更容易取得（比如图幅号对应的矩形）。

8.2 filters.crop：四种裁剪方式

filters.crop 支持 2D/3D 矩形、多边形、点+半径，以及它们的组合。选项：

8.2.1 bounds 的三种写法

bounds 是 PDAL 里最基础也最容易被写错的一个参数，它接受三种形式：

PDAL 经典写法
  3D：([xmin, xmax], [ymin, ymax], [zmin, zmax])
  2D：([xmin, xmax], [ymin, ymax])

GeoJSON BBOX
  3D：[xmin, ymin, zmin, xmax, ymax, zmax]
  2D：[xmin, ymin, xmax, ymax]

JSON 对象
  3D：{"minx":1,"miny":2,"minz":3,"maxx":101,"maxy":102,"maxz":103}
  2D：{"minx":1,"miny":2,"maxx":101,"maxy":102}

注意 PDAL 经典写法里每一维都是 [min, max]，而 GeoJSON BBOX 是“先全部 min 再全部 max”。这两种形式混着写不会报错，只会得到离谱的结果——比如把 [0,0,0,100,100,100] 当成 PDAL 写法解析，就变成 x∈[0,0]、y∈[0,100]、z∈[0,100] 的退化盒子，裁出 0 个点。

只有 JSON 对象形式支持额外带 crs 成员：

{
    "bounds": {
        "minx": 500000, "miny": 4000000,
        "maxx": 501000, "maxy": 4001000,
        "crs": "EPSG:32650"
    }
}

这条能力在跨坐标系裁剪时非常实用：数据是 EPSG:4326，但你手里的范围是 UTM 坐标，直接在 bounds 里写 crs 即可，不必先重投影。

8.2.2 多边形裁剪

[
    "input.las",
    {
        "type": "filters.crop",
        "polygon": "POLYGON((636889 851528, 636899 851475, 636976 851513, 636889 851528))"
    },
    "cropped.las"
]

polygon 可以给数组，每个多边形产出一个 PointView：

{
    "type": "filters.crop",
    "polygon": [
        "POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))",
        "POLYGON((200 200, 300 200, 300 300, 200 300, 200 200))"
    ]
}

配合 "output#.las" 就能一次裁出两份文件（第 5 章讲过这条组合）。

8.2.3 点 + 距离

{
    "type": "filters.crop",
    "point": "POINT(500000 4000000 100)",
    "distance": 500
}

二维点 POINT(x y) 裁出圆，三维点 POINT(x y z) 裁出球。做“局部样例抽取”很顺手：

pdal translate city.las sample.las -f crop \
    --filters.crop.point="POINT(121.47 31.23)" \
    --filters.crop.distance=500

注意 distance 单位与坐标单位一致。经纬度数据下 distance: 500 就是 500 度，会圈住整个地球——不报错，只是结果明显不对。

8.2.4 outside 与 a_srs

outside: true 反转裁剪逻辑，只保留范围外的点。典型用途是“抠掉一块”（比如挖去研究区内的水域）：

{
    "type": "filters.crop",
    "polygon": "POLYGON((500000 4000000, 501000 4000000, 501000 4001000, 500000 4001000, 500000 4000000))",
    "outside": true
}

a_srs 声明裁剪范围使用什么坐标系。不写时按点数据的 SRS 处理。两者不同时必须写，否则会静默裁错位置：

{
    "type": "filters.crop",
    "polygon": "POLYGON((120.1 30.1, 120.2 30.1, 120.2 30.2, 120.1 30.2, 120.1 30.1))",
    "a_srs": "EPSG:4326"
}

文档里有一条警告值得引用：如果同时存在多个不同 SRS 的 PointView，PDAL 会随机挑一个作为裁剪范围的 SRS 并打日志，这是不可靠的行为。多源输入时务必显式写 a_srs。

8.3 filters.reprojection：坐标变换

[
    {
        "filename": "input.las",
        "type": "readers.las",
        "spatialreference": "EPSG:26916"
    },
    {
        "type": "filters.reprojection",
        "in_srs": "EPSG:26916",
        "out_srs": "EPSG:4326"
    },
    "output.las"
]

选项表：

in_srs 是可省的，但建议总是写。省掉的后果是依赖数据自带的 SRS，而 LAS 文件里的 SRS 经常缺失或写错。显式指定相当于在 pipeline 里留了一条断言。

8.3.1 轴序：4326 的经典陷阱

EPSG:4326 的官方轴序是 纬度在前（lat, lon），而 GIS 数据习惯是 X=经度, Y=纬度。PROJ 6 之后更严格地遵守轴序定义，于是出现“重投影后 X/Y 互换”的现象。

PDAL 的处理方式：默认按常见习惯解释，需要时用轴序覆盖：

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

文档明确说明 "2, 1" 的作用是交换 X 与 Y，且明确提到“这对 EPSG:4326 这类坐标系常常是必需的”。所以：凡是输出到 4326，跑完必须验一次 X/Y 是否合理（经度应在 -180…180，纬度在 -90…90）。

8.3.2 重投影后必须重设 scale

这是第 2 章讲过的坑，在 pipeline 里更隐蔽。写 LAS 时坐标存成 int32，默认 scale = .01。UTM 坐标（如 500000）用 0.01 存没问题；重投影到 4326 后变成 116.39 这种值，0.01 的 scale 只能表达两位小数，横向误差以公里计。

正确写法（来自官方示例）：

[
    "input.las",
    {
        "type": "filters.reprojection",
        "in_srs": "EPSG:26916",
        "out_srs": "EPSG:4326"
    },
    {
        "type": "writers.las",
        "filename": "example-geog.las",
        "scale_x": "0.0000001",
        "scale_y": "0.0000001",
        "scale_z": "0.01",
        "offset_x": "auto",
        "offset_y": "auto",
        "offset_z": "auto"
    }
]

8.3.3 垂直基准与 geoidgrids

filters.reprojection 也能处理高程基准。官方示例是把椭球高（GRS80）转成正高（丹麦 DVR90）：

[
    "./1km_6135_632.laz",
    {
        "type": "filters.reprojection",
        "in_srs": "EPSG:25832",
        "out_srs": "+init=epsg:25832 +geoidgrids=C:/data/geoids/dvr90.gtx"
    },
    {
        "type": "writers.las",
        "filename": "output.laz",
        "a_srs": "EPSG:7416"
    }
]

要点：out_srs 用 PROJ 字符串挂 geoid 网格文件，写文件时再用 a_srs 把复合坐标系写进 VLR；不这么做，PDAL 会把垂直基准写成 “Unnamed Vertical Datum”。geoid 网格要事先备好，PROJ 不会自动下载。

8.4 filters.range 与 filters.expression：两种筛选语法

裁剪按几何筛点，这两个按属性筛点。它们功能重叠，但语法与未来走向不同。

8.4.1 filters.range：区间语法

[
    "input.las",
    {
        "type": "filters.range",
        "limits": "Z[0:100],Classification[2:2]"
    },
    "out.las"
]

区间写法规则（doc/stages/filters.range.md 的 Ranges 一节）：

[] 是闭区间（含端点），() 是开区间（不含端点），! 放在维度名后表示取反。

多条件的组合规则容易记错：

• 同一维度的多个区间 → OR
• 不同维度之间 → AND

文档给的例子最能说明：

Classification[1:2], Red[1:50], Blue[25:75], Red[75:255], Classification[6:7]

选中“分类为 1、2、6、7 之一，且蓝色在 25–75，且红色在 1–50 或 75–255”的点。分类被 OR 到一起，红蓝之间是 AND。

8.4.2 filters.expression：布尔表达式

[
    "input.las",
    {
        "type": "filters.expression",
        "expression": "0 <= Z <= 100 && Classification == 2"
    },
    "out.las"
]

表达式语法是 C 语言算子的子集，支持 *、/、+、-、!、>、>=、<、<=、==、!=、&&、||，运算顺序与 C 一致，支持括号。所有数学运算按 double 进行。

一条硬规则：没有自动的类型到布尔转换。下面这种写法是非法的：

((Intensity > 0) && Classification)      // 错
((Intensity > 0) && (Classification != 0))  // 对

官方文档在 filters.range 页面顶部给了一条建议：PDAL 2.5.x+ 起优先用 filters.expression，语法更简单，且与许多 stage 的 where 选项同源；filters.range 将在 PDAL 3.0 起废弃。

迁移对照：

8.4.3 where：筛选的第三种写法

第 5 章提过，每个 filter/writer 都支持 where。三者关系：

所以“只想让某些点参与这一步计算，但不删点”时，用 where；“这些点后面都不需要了”，才用 range/expression。

8.5 filters.assign：按条件改值

filters.assign 有两个选项，新写法用 value：

[
    "input.las",
    {
        "type": "filters.assign",
        "value": [
            "Classification = 2 WHERE HeightAboveGround < 5",
            "Classification = 1 WHERE HeightAboveGround >= 5"
        ]
    },
    "out.laz"
]

语法是 维度 = 值表达式 [WHERE 条件表达式]。值表达式支持 +、-、*、/ 与常量、维度名。

常用场景是颜色位深换算：Red = Red / 256 把 16 位降到 8 位；反过来 Red = Red * 256 是给 COPC（LAS 1.4 要求 16 位色）准备数据，三个通道都要写进 value 数组。

两个必须知道的点：

1. assignment 与 condition 已废弃，虽然还能用。新代码用 value。
2. PDAL 2.7.0 起，给不存在的维度赋值会自动创建它，类型固定 double。 这条变化让“先 ferry 再 assign”的老写法不再必需，但要注意自动创建出来的是 double，写进 LAS 时若声明成 uint16 会有一次转换。

8.6 filters.ferry：在维度之间搬数据

filters.ferry 的选项只有一个：dimensions，格式 <from>=><to>，多个用逗号分隔，空格被忽略。

{
    "type": "filters.ferry",
    "dimensions": "X => StatePlaneX, Y => StatePlaneY"
}

两个用法：

保存重投影前的坐标——重投影会就地改写 X/Y，想留底就得先 ferry：

[
    "input.las",
    { "type": "filters.ferry", "dimensions": "X => OrigX, Y => OrigY" },
    { "type": "filters.reprojection", "out_srs": "EPSG:4326" },
    "out.las"
]

创建新维度——from 留空：

[
    { "type": "readers.gdal", "filename": "somefile.tif" },
    { "type": "filters.ferry", "dimensions": "=>Classification" },
    { "type": "filters.assign", "assignment": "Classification[:]=2" },
    "out.las"
]

注意例子里的 assignment 是旧语法（Classification[:]=2），新写法应为 "value": "Classification = 2"。

to 维度不存在时会自动创建，初值 0；from 必须存在，否则报错。旧语法用 = 代替 => 也仍然支持。

8.7 filters.transformation：任意 4×4 变换

当需要的变换不是“换坐标系”而是“挪一下、转一下、缩放一下”时，用 filters.transformation：

{
    "type": "filters.transformation",
    "matrix": "0 -1 0 1  1 0 0 2  0 0 1 3  0 0 0 1"
}

matrix 是 16 个空格分隔的数，行主序。四种基本形式：

平移 (tx, ty, tz)
1 0 0 tx   0 1 0 ty   0 0 1 tz   0 0 0 1

绕 Z 轴旋转 θ
cosθ -sinθ 0 0   sinθ cosθ 0 0   0 0 1 0   0 0 0 1

缩放 (sx, sy, sz)
sx 0 0 0   0 sy 0 0   0 0 sz 0   0 0 0 1

两个要点：

• 它完全不考虑坐标系。文档明确写了，这个 filter 不应用也不考虑任何 SRS 信息。用它做“坐标系纠正”是错的，那是 reprojection 的活。
• invert: true 会应用矩阵的逆，适合“先按 A 变换、之后要还原”的场景。

8.8 filters.overlay：用多边形给点赋属性

裁剪是“留不留”，overlay 是“改成什么”。它读一个 OGR 数据源，把落在多边形内的点赋上该多边形的属性值：

[
    "autzen-dd.las",
    {
        "type": "filters.overlay",
        "dimension": "Classification",
        "datasource": "attributes.shp",
        "layer": "attributes",
        "column": "CLS"
    },
    "attributed.las"
]

典型用途：把土地利用矢量里的类别码写成点的 Classification，或把行政区代码写成自定义维度。注意它只做空间包含判断，不处理坐标系差异——数据源与点云的 SRS 必须一致，否则要先用 ogr 查询里的 t_srs 或提前转换。

8.9 组合实例：一条完整的裁剪 + 重投影管道

把本章的 stage 串起来，做一个常见任务：从大范围 LAZ 里取一块研究区，统一到 UTM，并按地面/非地面分两份输出。

[
    {
        "type": "readers.las",
        "filename": "province.laz",
        "override_srs": "EPSG:4490"
    },
    {
        "type": "filters.expression",
        "expression": "Z >= 0 && Z <= 3000"
    },
    {
        "type": "filters.crop",
        "polygon": "POLYGON((120.10 30.10, 120.30 30.10, 120.30 30.30, 120.10 30.30, 120.10 30.10))",
        "a_srs": "EPSG:4326"
    },
    {
        "type": "filters.reprojection",
        "in_srs": "EPSG:4490",
        "out_srs": "EPSG:32650"
    },
    {
        "type": "writers.las",
        "filename": "study_area.laz",
        "compression": "true",
        "scale_x": 0.001,
        "scale_y": 0.001,
        "scale_z": 0.001,
        "offset_x": "auto",
        "offset_y": "auto",
        "offset_z": "auto",
        "a_srs": "EPSG:32650"
    }
]

顺序上的三个判断：

1. expression 放在 crop 之前。Z 的异常值先在原坐标筛掉，比裁完再筛省一次比对；
2. crop 用 a_srs 声明自己的坐标系，因为数据是 4490（地理坐标）、裁剪范围也是 4326 经纬度，两者轴序与基准都可能有差异，显式声明是唯一可靠做法；
3. reprojection 放在 crop 之后。少变换一批点。

readers.las 上的 override_srs 是第一道保险：文件头 SRS 不可信时用它钉死；省掉时要在日志里确认 SRS 被正确识别。

8.10 排错

error_on_failure 那条值得展开：默认 false 意味着投影失败的点会被静默跳过（保留原值或置零），批量处理时你会得到一份“看起来正常、局部有洞”的结果。生产任务里把它设成 true，让失败显式暴露，比事后排查便宜得多。

8.11 生产建议

先裁后投。 除非裁剪范围只在目标坐标系下才方便表达，否则裁剪一定放在重投影之前。省的不只是重投影的计算，还有中间数据的体积。

坐标范围一律显式声明 SRS。 filters.crop.a_srs、readers.*.override_srs、writers.*.a_srs 三处各管一段。三个都写清楚了，跨系统交付时就不会出现“在我这看着对、到他那偏了 500 米”。

重投影后立刻验坐标范围。 用 pdal info --summary 看 minx/maxx 是否落在预期区间（4326 应落在 ±180/±90 内）。这一步能挡住轴序与单位两类错误。

把 scale 与 SRS 绑定管理。 地理坐标必须配 1e-7 级 scale，这条规则写进模板，别靠人记。

用 filters.expression 替代 filters.range。 后者将在 3.0 废弃，新 pipeline 直接用表达式，语法与 where 一致，迁移成本一次性付清。

8.12 小结