用LiDAR点云估算楼高

LiDAR 数据可以非常强大地提取城市规模的地面和物体的高程。在 One Concern，我们正在使用 LiDAR 数据提取地面和建筑物高程，以改进进入我们的自然灾害模型的暴露信息，最终估计洪水和地震的影响。

由于 3DEP 项目预计将在 2023 年之前收集全美 LiDAR 数据，因此 LiDAR 数据的可用性应该会变得更加容易和广泛得多。然而，由于 LiDAR 是点云数据，每个城市可能包含数十亿个点，读取和处理数据并不简单。尽管已经有其他文章介绍如何处理 LiDAR 数据以提取建筑物高度，但其中一些库已经停用或被取代。在这篇文章中，我打算提供一组更新的工具和流程，以使用 LiDAR 数据提取建筑物高程，以旧金山为例。

1、工具

我将使用以下工具——

pdal — 用于读取 LiDAR 文件
带有点云的 PostgreSQL — 用于存储 LiDAR 点云数据的数据库和扩展
docker — 在虚拟环境中设置工具（可选）
QGIS 3.x — 可视化数据（可选）

2、数据

对于本文，我将使用两个数据集——LiDAR 和 OSM。

2.1 激光雷达数据

旧金山地区的 LiDAR 数据可从 USGS 数据库获得，在该数据库中它被称为 Golden Gate LiDAR 项目 . 该数据是 2010 年为旧金山州立大学收集的，覆盖了旧金山及其附近的 835 平方英里。

你可以通过 LASzip 网站上编译的链接访问其他一些公开可用的全球 LiDAR 数据集。

2.2 OpenStreetMaps 数据

为了将 LiDAR 数据分配给建筑足迹，我将使用来自 OpenStreetMaps (OSM) 的北加州建筑足迹形状文件 (.shp.zip) 数据。

3、设置环境

我使用 docker 在虚拟容器中设置我的环境。我在下面使用 docker 为 PDAL 和 PostgreSQL 提供了使用 pointcloud 的设置步骤。虽然我在 MacOS 上使用 docker，但它也适用于 Linux 和 Windows，并且无论你的操作系统如何，设置都应该保持不变。你还可以使用各自网站上的说明直接在你的操作系统上使用 pointcloud 设置 PDAL 和 PostgreSQL。

如果你使用的是 docker，最简单的方法是首先为容器创建一个桥接网络以相互通信。以下命令创建一个名为 pdal-net 的桥接网络，但是你可以将其替换为你想要的名称。

docker network create pdal-net

我将使用 PDAL 读取我的 LiDAR 文件。 PDAL 是一个类似 GDAL 的库，旨在用于点云数据。它是以前用于相同目的的库（如 libLAS 和 LASlib）的新替代品。在您的终端中，从包含 LiDAR 数据的当前目录，使用以下命令在后台启动名为 pdal 的 pdal docker 容器，并将其附加到网络 pdal-net。标志 -v $(pwd):/scripts 在 pdal 容器的 /scripts 目录中的主机系统上安装你的当前目录。

docker run --network pdal-net -v $(pwd):/scripts --name pdal pdal/pdal:latest

我将使用带有点云扩展的 PostgreSQL 来存储和处理 LiDAR 数据。由于数据包含数百万个点，将其存储在数据库中更易于处理，点云扩展提供了一种使用 PostgreSQL 查询存储、访问和处理数据的便捷方式。在与上述相同的目录中，使用以下命令在名为 pgpointcloud 的 docker 容器中设置并启动数据库，其中包含必要的 postgis、pointcloud 和 pointcloud_postgis 扩展。

docker run --network pdal-net --name pgpointcloud -e POSTGRES_DB=pointclouds -e POSTGRES_PASSWORD=mysecretpassword -v $(pwd):/mnt -d pgpointcloud/pointcloud:latest

4、读取激光雷达数据

LiDAR 数据通常以 .las 或 .laz 文件的形式提供，其中后者是占用较少磁盘空间的压缩版本。旧金山数据以 .las 和 .laz 文件形式提供。在 PDAL 中读取这两种类型文件的步骤是相同的，只要安装了正确的 LiDAR 压缩库依赖项（如 LASzip），如果你使用上述 docker 镜像，它们就是这样。 LiDAR 数据通常收集在统一的网格中，每个网格的数据作为单独的 .las/.laz 文件提供。旧金山的数据被分成近 1200 个网格，压缩文件占用近 38GB 的磁盘空间。

在 PDAL 中，读取、过滤和写入点云文件等处理步骤在称为管道的 JSON 脚本中进行描述。以下管道读取其中一个 LiDAR 文件并将它们输出为 .csv（逗号分隔文本）文件，你可以将其导入 QGIS 以进行可视化。

{
  "pipeline":[
    {
      "type":"readers.las",
      "filename":"laz/ARRA-CA_GoldenGate_2010_001019.laz",
      "spatialreference":"EPSG:26910"
    },
    {
      "type":"writers.text",
      "format":"csv",
      "keep_unspecified":"true",
      "filename":"laz/ARRA-CA_GoldenGate_2010_001019.laz.txt"
    }
  ]
}

上述代码使用las reader读取laz文件夹下的ARRA-CA_GoldenGate_2010_001019.laz文件，使用text writer输出对应的.csv文件。根据此数据集提供的元数据，我确定该数据是在空间参考 EPSG:26910 中提供的，然后将其指定给读者。

要运行上面的管道代码，请将上面的文本作为 pipeline_single.json 保存在包含 laz 文件夹的目录中。然后我们将从 pdal 容器中运行管道。使用 Bash 进入 pdal 容器
docker 启动 pdal && docker exec -it pdal bash。接下来，转到容器中包含 pipeline_single.json 文件 (cd /scripts) 的正确目录，最后运行以下命令，这将在 laz 文件夹中创建所需的 .csv 文件。

pdal pipeline -i pipeline_single.json

.csv 文件中的文本格式的点云数据可以导入 QGIS 并在下面可视化。

虽然文本文件可以更轻松地可视化 LiDAR 数据，但它们并不是存储和处理整个城市的大量数据的理想媒介。我们将在 PostgreSQL 中做到这一点。 PostgreSQL 提供数据库功能，postgis 扩展使处理几何和地理数据变得容易，点云扩展提供了一种存储和提取点云数据的便捷方式。 PDAL 有一个处理点云数据的 postgres 编写器，因此我们可以使用它来读取 .laz 文件并将它们直接存储到我们的数据库中。

如果你如上所述启动了 pgpointcloud 的 docker 容器，那么 PostgreSQL 数据库应该已经运行了相关的扩展。你可以使用以下 PDAL 管道代码读取数据集的所有 LiDAR 数据文件并将它们存储在您的数据库中。将以下 PDAL 管道代码保存在名为 pipeline_batch.json 的文件中。

{
  "pipeline":[
    {
      "type":"readers.las",
      "spatialreference":"EPSG:26910"
    },
    {
      "type":"filters.chipper",
      "capacity":600
    },
    {
      "type":"writers.pgpointcloud",
      "connection":"host='pgpointcloud' dbname='pointclouds' user='postgres' password='mysecretpassword' port='5432'",
      "compression":"dimensional",
      "srid":"26910"
    }
  ]
}

请注意，在上面的代码中，LiDAR 数据的文件名未在管道中指定，因为我们希望管道能够从目录中读取不同的文件。调用管道时，文件名将通过 bash 脚本传递给管道。此外，还有一个容量为 600 的削片（chipper）过滤器。它将点云数据组合成 600 个附近点的补丁（patch），数据库存储这些补丁而不是单个点。我将在下面进一步讨论补丁的概念。建议单个补丁中包含 400 到 600 个点，以使每个补丁保持在 8kb 的 PostgreSQL 页面大小限制之下。最后，此管道中的编写器连接到 PostgreSQL 数据库，而不是像在之前的管道中那样写入 .csv 文件。

以下 bash 脚本查找指定目录中的所有 .laz 文件，并将每个文件作为参数传递给 PDAL 管道，然后 PDAL 管道读取文件并将其存储在名为 lidar 的 PostgreSQL 表中。在我的实施中，我首先只选择了与旧金山收集网格相关的文件子集，并将它们存储在子目录 san francisco 中，以减少要读取和存储的文件数量。

#!/bin/bash
dirname="/scripts/LIDAR/ARRA-CA_GoldenGate_2010"
subdirname="san francisco"
tablename="lidar"
for filename in "$dirname"/"$subdirname"/*.laz; do
  echo "Processing $filename ..."
  pdal pipeline -i "$dirname"/pipeline_batch.json \
  --readers.las.filename="$filename" \
  --writers.pgpointcloud.table="$tablename"
done

将上述文件另存为 batchLAZtoPostgres.sh，使其可执行，并使用以下命令从 pdal 容器中运行它。请注意，/scripts/LIDAR/ARRA-CA_GoldenGate_2010 是 pdal 容器中 LiDAR 文件的目录结构。

# make script executable
chmod +x batchLAZtoPostgres.sh
# run the script
./batchLAZtoPostgres.sh

由于 pdal 和 pgpointcloud 容器都在同一个桥接网络上，PDAL 将能够连接到数据库并将点云数据存储为补丁。根据文件的数量和大小，读取和存储数据库中的所有点可能需要几小时到几天的时间。

5、在 PostgreSQL 中处理 LiDAR 数据

补丁是一种在 PostgreSQL 表中存储点云数据的有效方法，其中每个补丁都包含许多附近的点。这有助于大大减少数据的行数，并使对表的操作更快。补丁和点的模式在数据库中的表 pointcloud_formats 中描述。此外，点云扩展提供了同时处理点和面片的功能，包括为每个点或面片的边界框范围提取单独的属性，如高程、强度或坐标。如果想阅读更多内容，OpenGEO 的 Paul Ramsey 的这篇演讲对点和补丁的示例提供了非常好的解释。

所有后续步骤都在 pgpointcloud 容器中实现，你可以使用以下命令对其进行 bash -

docker start pgpointcloud
docker exec -it pgpointcloud bash

我们将把 LiDAR 数据与 OSM 建筑物足迹结合起来，以提取建筑物的高程和高度。因此，让我们首先从 pgpointcloud 容器中使用以下命令将上述 OSM shapefile 数据加载到名为 osm 的表中 -

shp2pgsql -s 4326 gis_osm_buildings_a_free_1.shp public.osm | psql -U postgres -d pointclouds

我们将在其余大部分步骤中使用 PostgreSQL 查询。可以使用以下命令访问数据库并从数据库中运行所有查询。

psql -U postgres -d pointclouds

现在，我们准备提取建筑物的立面和高度，详细过程如下所述。获取建筑物高程的概览过程是首先从 OSM 识别建筑物覆盖区内的所有 LiDAR 点，然后从这些点集合中获取建筑物高程。获取建筑物高度的概述是首先沿着每个建筑物占地面积的轮廓估计地面高程，然后从建筑物高程中减去它。

5.1 创建Lidar空间索引

预处理激光雷达表以创建空间索引

a. 查看补丁中每个点的属性。补丁和点的格式在名为 pointcloud_formats 的单独表中以 XML 格式描述。

SELECT * FROM pointcloud_formats;

b. 查看表中的列表。您会注意到补丁存储在 pa 列中。

\d lidar

c. 查看表中的补丁数和总点数——

SELECT Count(*), Sum(PC_NumPoints(pa)) FROM lidar;

d. 在激光雷达表上创建空间索引以加快查询速度（观察到的加速为 40 倍）—

CREATE INDEX lidar_envelope_pkey ON lidar USING GIST(PC_EnvelopeGeometry(pa));

5.2 创建OSM空间索引

预处理 osm 表以创建空间索引

a. 可以使用下面的查询删除城市边界框外的点。将 (lng_min, lat_min, lng_max, lat_max) 替换为感兴趣的边界框的范围，其中 lng 指经度，lat 指纬度。

DELETE FROM osm WHERE osm.gid IN (SELECT a.gid FROM osm a WHERE NOT ST_Intersects(a.geom, ST_MakeEnvelope (lng_min, lat_min, lng_max, lat_max, 4326)) );

b. 我建议将 osm geomtery 列 geom 转换为与激光雷达表相同的 SRID（空间参考系统）。 UTM 区域非常适合处理小区域和测量，因此我们将 geom 转换为 EPSG:26910，这与 LiDAR 数据的参考相同。

ALTER TABLE osm ALTER COLUMN geom TYPE geometry(MultiPolygon, 26910) USING ST_Transform(geom, 26910);
CREATE INDEX osm_26910_pkey ON osm USING GIST (geom);

5.3 创建点补丁列

在 osm 中创建新列以包含建筑物足迹内的点补丁

a. 为了向 osm 添加一个新列来存储补丁，首先使用 \d lidar 检查激光雷达中的补丁类型。通常是 pcpatch(1)。在 osm 中添加一个新的列来存储补丁

ALTER TABLE osm ADD COLUMN pa pcpatch(1);

b. 存储覆盖每个建筑物足迹多边形的点云补丁。下面的查询首先识别与每个建筑物足迹相交的所有补丁，然后计算它们的并集，最后找到并集补丁与足迹的交集。预计此查询需要几个小时。

UPDATE osm SET pa = sq.pa FROM (WITH patches AS (SELECT o.gid AS gid, o.geom AS geom, l.pa AS pa FROM lidar AS l JOIN osm AS o ON PC_INTERSECTS(l.pa, o.geom)) SELECT gid, PC_INTERSECTION(PC_UNION(pa), geom) AS pa FROM patches GROUP BY gid, geom) AS sq WHERE osm.gid = sq.gid;

5.4 计算统计信息

我们将建筑物标高描述为建筑物占地面积内所有点的标高的统计数据。可以根据你的用例选择统计信息。在这里，我计算了所有点的最大、平均和中值海拔。请注意，LiDAR 数据中的 Z 值是给定元数据中指定的特定数据（在本例中为 NAD83）的高程。不是直接取最大值，而是将其估计为第 99.9 个百分位值，以减少异常值导致高估最高海拔的可能性。

a. 将相关列添加到 osm —

ALTER TABLE osm ADD COLUMN z_avg DOUBLE PRECISION NULL, ADD COLUMN z_median DOUBLE PRECISION NULL, ADD COLUMN z_max DOUBLE PRECISION NULL;

b. 计算和存储高程统计数据

UPDATE osm SET z_avg = sq.z_avg, z_median = sq.z_median, z_max = sq.z_max FROM (WITH patches AS (SELECT o.gid AS gid, PC_GET(PC_EXPLODE(o.pa), ‘Z’) AS pt_z FROM osm AS o) SELECT gid, AVG(pt_z) AS z_avg, PERCENTILE_CONT(0.5) WITHIN GROUP(ORDER BY pt_z) AS z_median, PERCENTILE_CONT(0.999) WITHIN GROUP(ORDER BY pt_z) AS z_max FROM patches GROUP BY gid) AS sq WHERE osm.gid = sq.gid;

5.5 获取地面高程

通过识别建筑物轮廓并将其与 LiDAR 数据相结合来获取地面高程。

a. 创建一个名为 osm_outline 的新表作为 osm 表的副本，以处理通过获取建筑物多边形周围 2m 和 1m 缓冲区之间的差异生成的建筑物轮廓。

CREATE TABLE osm_outline AS SELECT gid, osm_id, geom FROM osm;
UPDATE osm_outline SET geom = buffer FROM (SELECT o.gid, ST_MULTI(ST_DIFFERENCE(ST_MULTI(ST_Buffer(o.geom, 2)), ST_MULTI(ST_BUFFER(o.geom, 1)))) AS buffer FROM osm_outline AS o) AS sq WHERE osm_outline.gid = sq.gid;
CREATE INDEX osm_outline_pkey ON osm_outline USING GIST (geom);

b. 存储与轮廓相交的点云补丁。下面的查询首先识别与每个轮廓相交的所有补丁，然后计算它们的并集，最后找到并集补丁与轮廓的交集。预计此查询需要几个小时。

ALTER TABLE osm_outline ADD COLUMN pa pcpatch(1);
UPDATE osm_outline SET pa = sq.pa FROM (WITH patches AS (SELECT o.gid AS gid, o.geom AS geom, l.pa AS pa FROM lidar AS l JOIN osm_outline AS o ON PC_INTERSECTS(l.pa, o.geom)) SELECT gid, PC_INTERSECTION(PC_UNION(pa), geom) AS pa FROM patches GROUP BY gid, geom) AS sq WHERE osm_outline.gid = sq.gid;

c. 存储与建筑物轮廓相交的点云的最小 Z 值。不是直接取最小值，而是将其估计为第一个百分位值，以减少异常值导致低估地面高程的可能性。

ALTER TABLE osm_outline ADD COLUMN z_ground DOUBLE PRECISION NULL;
UPDATE osm_outline SET z_ground = sq.z_min FROM (WITH patches AS (SELECT o.gid AS gid, PC_GET(PC_EXPLODE(o.pa), ‘Z’) AS pt_z FROM osm_outline AS o) SELECT gid, PERCENTILE_CONT(0.01) WITHIN GROUP(ORDER BY pt_z) AS z_min FROM patches GROUP BY gid) AS sq WHERE osm_outline.gid = sq.gid;

d. 将地面高程添加到原始 osm 表，并将高度计算为脚印高程和地面高程之间的差值。

ALTER TABLE osm ADD COLUMN z_ground DOUBLE PRECISION NULL, ADD COLUMN height_avg DOUBLE PRECISION NULL, ADD COLUMN height_median DOUBLE PRECISION NULL, ADD COLUMN height_max DOUBLE PRECISION NULL;
UPDATE osm AS o SET z_ground = oo.z_ground FROM osm_outline oo WHERE o.gid = oo.gid;
UPDATE osm SET height_avg = (z_avg — z_ground), height_median = (z_median — z_ground), height_max = (z_max — z_ground);

就是这样！我们成功计算出每个建筑物的建筑物标高、地面标高和建筑物高度，并且所有这些数据连同每个建筑物的建筑物占地面积几何图形一起存储在 osm 表中。

6、可视化

如果想可视化海拔高度，你可以将数据导出到 Shapefile。首先使用 \q 退出数据库，然后运行以下命令将 osm 表导出到名为 sanfrancisco_buildings.shp 的文件 —

pgsql2shp -f sanfrancisco_buildings.shp -g geom -u postgres -d pointclouds “SELECT gid, osm_id, code, fclass, name, type, geom, z_avg, z_median, z_max, z_ground, height_avg, height_median, height_max FROM osm WHERE height_median IS NOT NULL”

我将 Shaepefile 导入 QGIS，并使用 Qgis2threejs 插件创建交互式输出，如本文顶部所示。该插件允许以 3D 形式显示地面高程和建筑物，并且可以从任何浏览器轻松访问可视化。

原文链接：Estimating Building Heights Using LiDAR Data

BimAnt翻译整理，转载请标明出处