pdal 3.4.5

PyPI

PDAL Python支持可通过PyPI安装

pip install PDAL

GitHub

PDAL Python扩展的仓库可在https://github.com/PDAL/python找到

从PDAL 1.7版本开始，Python支持独立于PDAL本身发布。

简单

以下是一个简单的管道示例，该管道简单地读取一个ASPRS LAS文件并按X维度对其进行排序

json = """
{
  "pipeline": [
    "1.2-with-color.las",
    {
        "type": "filters.sort",
        "dimension": "X"
    }
  ]
}"""

import pdal
pipeline = pdal.Pipeline(json)
count = pipeline.execute()
arrays = pipeline.arrays
metadata = pipeline.metadata
log = pipeline.log

程序化管道构建

上一个示例指定了管道为JSON字符串。另外，可以通过创建Stage实例并将它们管道连接起来来构建管道。例如，上一个管道可以指定为

pipeline = pdal.Reader("1.2-with-color.las") | pdal.Filter.sort(dimension="X")

>>> help(pdal.Filter.head)
Help on function head in module pdal.pipeline:

head(**kwargs)
    Return N points from beginning of the point cloud.

    user_data: User JSON
    log: Debug output filename
    option_file: File from which to read additional options
    where: Expression describing points to be passed to this filter
    where_merge='auto': If 'where' option is set, describes how skipped points should be merged with kept points in standard mode.
    count='10': Number of points to return from beginning.  If 'invert' is true, number of points to drop from the beginning.
    invert='false': If true, 'count' specifies the number of points to skip from the beginning.

# update pipeline in-place
pipeline = pdal.Pipeline()
pipeline |= stage
pipeline |= pipeline2

使用Numpy数组读取

以下更复杂的场景演示了PDAL和Python之间完整的循环

• 从GitHub读取一个小测试文件到Numpy数组

• 使用Numpy过滤数组以获取强度

• 将过滤后的数组传递给PDAL再次进行过滤

• 使用TileDB-PDAL集成和TileDB writer插件将最终的过滤数组写入LAS文件和一个TileDB数组

import pdal

data = "https://github.com/PDAL/PDAL/blob/master/test/data/las/1.2-with-color.las?raw=true"

pipeline = pdal.Reader.las(filename=data).pipeline()
print(pipeline.execute())  # 1065 points

# Get the data from the first array
# [array([(637012.24, 849028.31, 431.66, 143, 1,
# 1, 1, 0, 1,  -9., 132, 7326, 245380.78254963,  68,  77,  88),
# dtype=[('X', '<f8'), ('Y', '<f8'), ('Z', '<f8'), ('Intensity', '<u2'),
# ('ReturnNumber', 'u1'), ('NumberOfReturns', 'u1'), ('ScanDirectionFlag', 'u1'),
# ('EdgeOfFlightLine', 'u1'), ('Classification', 'u1'), ('ScanAngleRank', '<f4'),
# ('UserData', 'u1'), ('PointSourceId', '<u2'),
# ('GpsTime', '<f8'), ('Red', '<u2'), ('Green', '<u2'), ('Blue', '<u2')])
arr = pipeline.arrays[0]

# Filter out entries that have intensity < 50
intensity = arr[arr["Intensity"] > 30]
print(len(intensity))  # 704 points

# Now use pdal to clamp points that have intensity 100 <= v < 300
pipeline = pdal.Filter.range(limits="Intensity[100:300)").pipeline(intensity)
print(pipeline.execute())  # 387 points
clamped = pipeline.arrays[0]

# Write our intensity data to a LAS file and a TileDB array. For TileDB it is
# recommended to use Hilbert ordering by default with geospatial point cloud data,
# which requires specifying a domain extent. This can be determined automatically
# from a stats filter that computes statistics about each dimension (min, max, etc.).
pipeline = pdal.Writer.las(
    filename="clamped.las",
    offset_x="auto",
    offset_y="auto",
    offset_z="auto",
    scale_x=0.01,
    scale_y=0.01,
    scale_z=0.01,
).pipeline(clamped)
pipeline |= pdal.Filter.stats() | pdal.Writer.tiledb(array_name="clamped")
print(pipeline.execute())  # 387 points

# Dump the TileDB array schema
import tiledb
with tiledb.open("clamped") as a:
    print(a.schema)

执行可流式传输的管道

可流式传输的管道（仅由可流式传输的PDAL阶段组成的管道）可以通过Pipeline.iterator()以流式模式执行。这返回一个迭代器对象，每次产生最多chunk_size大小（默认=10000）的Numpy数组。

import pdal
pipeline = pdal.Reader("test/data/autzen-utm.las") | pdal.Filter.range(limits="Intensity[80:120)")
for array in pipeline.iterator(chunk_size=500):
    print(len(array))
# or to concatenate all arrays into one
# full_array = np.concatenate(list(pipeline))

Pipeline.iterator() 还可以接受一个可选的 prefetch 参数（默认=0），允许并行预取多达这个数量的数组，并在返回给调用者之前进行缓冲。

如果您只想以流式模式执行可流式处理的流水线，并且不需要访问数据点（通常流水线包含 Writer 阶段时），可以使用 Pipeline.execute_streaming(chunk_size) 方法。这与 sum(map(len, pipeline.iterator(chunk_size))) 功能上等效，但更高效，因为它避免了在内存中分配和填充任何数组。

访问网格数据

一些 PDAL 阶段（例如 filters.delaunay）创建 TIN 类型的网格数据。

可以使用 Pipeline.meshes 属性在 Python 中访问这些数据，该属性返回一个形状为 (1,n) 的 numpy.ndarray，其中 n 是网格中的三角形数量。

如果 PointView 中不包含网格数据，则 n = 0。

每个三角形是一个元组 (A,B,C)，其中 A、B 和 C 是指向 PointView 的索引，用于标识构成三角形的顶点。

Meshio 集成

meshes 属性提供了面数据，但不易作为网格使用。因此，我们提供了可选的 Meshio 库集成。

pdal.Pipeline 类提供了一个 get_meshio(idx: int) -> meshio.Mesh 方法。该方法从 PointView 数组和网格属性创建一个 Mesh 对象。

功能简单使用示例

import pdal

...
pl = pdal.Pipeline(pipeline)
pl.execute()

mesh = pl.get_meshio(0)
mesh.write('test.obj')

高级网格使用案例

USE-CASE：取一个激光雷达地图，从地面点创建网格，分割成瓦片，并将瓦片存储在 PostGIS 中。

(示例使用 1.2-with-color.las 并未进行地面分类以便清晰)

import pdal
import psycopg2
import io

pl = (
    pdal.Reader(".../python/test/data/1.2-with-color.las")
    | pdal.Filter.splitter(length=1000)
    | pdal.Filter.delaunay()
)
pl.execute()

conn = psycopg(%CONNNECTION_STRING%)
buffer = io.StringIO

for idx in range(len(pl.meshes)):
    m =  pl.get_meshio(idx)
    if m:
        m.write(buffer,  file_format = "wkt")
        with conn.cursor() as curr:
          curr.execute(
              "INSERT INTO %table-name% (mesh) VALUES (ST_GeomFromEWKT(%(ewkt)s)",
              { "ewkt": buffer.getvalue()}
          )

conn.commit()
conn.close()
buffer.close()