PDAL в Python: Практическое руководство по конвейерам обработки облаков точек

Узнайте, как запускать конвейеры PDAL в Python, обрабатывать облака точек LAZ, обрезать и объединять данные, а также создавать GeoTIFF-файлы ЦМР из точек, классифицированных как поверхность грунта.

PDAL (Point Data Abstraction Library) — это открытая библиотека и набор инструментов для обработки облаков точек. Она считается своего рода “GDAL для облаков точек” — так же, как GDAL работает с растровыми и векторными геоданными, PDAL работает с облаками точек. Я предпочитаю использовать PDAL для решения всех задач обработки облаков точек. PDAL прекрасно работает и с данными LiDAR и с облаками точек получеными методом фотограмметрии.

В этой статье рассмотрим несколько полезных и простых примеров использования PDAL, которые позволят вам использовать вам PDAL для решения ваших непосредственных задач, или интегрировать в ваши пайплайны обработки данных. Так же расмотрим пример создания цифровой модели рельефа в формате GeoTIFF на основании LAZ файла. В тексте статьи приведены блоки кода которые вы можете просто копировать и адаптировать под ваши задачи.

Навигация по статье

• Пайплайны
• Базовый пример запуска PDAL пайплайна
• Запуск PDAL пайплайна в Python
• Чтение Облака точек в NumPy array
• Простые но полезные пайплайны
• Обрезка облака точек по экстенту
• Обрезка облака точек по полигону
• Merging нескольких point clouds
• Создание GeoTIFF DEM из облака точек
• Подведение итогов

Пайплайны

Главный принцип работы PDAL, который может напугать новичка — это pipeline (конвейер обработки данных), который описывается в виде JSON-структуры.

JSON пайплайн представляет из себя список в котором последовательно перечислены следующие элементы:

Input (Reader) → Filters → Output (Writer)

Компоненты pipeline:

1. Readers (чтение)

   Загружают данные (например, LAS/LAZ файлы)

2. Filters (обработка)

   Преобразуют данные:

   • фильтрация шума
   • классификация
   • merging
   • clipping
   • кластеризация
   • и много других инструментов filters.
3. Writers (запись)

   Сохраняют результат (в файл или базу данных)

Каждый шаг — независимый блок:

• легко можно менять порядок
• можно их комбинировать

Базовый пример запуска PDAL пайплайна

Простейший способ это записать пайплайна в отдельный JSON-файл и потом использовать

Записываем наш пайплайн в отдельный JSON-файл my_pipeline.json

{
  "pipeline": [
    "input.laz",
    {
      "type": "filters.outlier",
      "method": "statistical",
      "mean_k": 8,
      "multiplier": 2.0
    },
    {
      "type": "filters.smrf"
    },
    {
      "type": "filters.expression",
      "expression": "Classification == 2"
    },
    "output.laz"
  ]
}

Далее в терминале выполняем команду:

pdal pipeline my_pipeline.json

Что здесь происходит:

• читаем файл input.laz
• удаляем выбросы (шум)
• классифицируем землю при помощи Simple Morphological Filter (SMRF)
• оставляет только ground points (Classification == 2)
• сохраняем результат в output.laz

получаем такой результат:

Запуск PDAL пайплайна в Python

Рассмотрим этот же пайплайн но запущенный в Python и пути к файлам передадим через переменные.

import json  # нужен, чтобы превратить Python-структуру данных в JSON-строку
import pdal  # Python binding для PDAL, через него создаётся и запускается pipeline.

# Пути к файлам
input_laz_file = "input.laz"
output_laz_file = "output.laz"

# Описание pipeline
my_pipeline = [
    # Этап 1 — чтение входного файла
    input_laz_file,
    # Этап 2 — удаление выбросов
    {
        "type": "filters.outlier",  # Фильтр для поиска outliers — подозрительных точек, слишком удалённых от соседей.
        "method": "statistical",  # Используется статистический метод
        "mean_k": 8,  # Для каждой точки рассматриваются 8 соседей
        "multiplier": 2.0  # Порог “насколько далеко от нормы” считать точку выбросом.
    },
    # Этап 3 — классификация земли через SMRF
    {
        "type": "filters.smrf"
    },
    # Этап 4 — оставить только ground points
    {
        "type": "filters.expression",
        "expression": "Classification == 2"
    },
    # Этап 5 — запись результата
    {
        "type": "writers.las",
        "filename": output_laz_file
    }
]

# Преобразование pipeline в JSON
pipeline_json = json.dumps(my_pipeline)

# Просто отладочный вывод. Можем распечатать и посмотреть коректно все ли коректно.
# Полезно, чтобы:
# - проверить итоговую структуру pipeline
# - убедиться, что JSON собрался корректно
# - при необходимости скопировать его и запустить отдельно через CLI
print(f"pipeline_json: {pipeline_json}")

# Здесь создаётся Python-объект PDAL pipeline.
# На этом этапе pipeline ещё не выполняется.
pipeline = pdal.Pipeline(pipeline_json)

# Выполнение pipeline
# возвращает число обработанных точек
count = pipeline.execute()

# метаданные выполнения pipeline
metadata = pipeline.metadata

# текстовый лог выполнения. Если пайплайн выполнился благополучно, то log пустой.
log = pipeline.log

print(f"metadata: {metadata}")
print(f"count: {count}")
print(f"log: {log}")

Чтение Облака точек в NumPy array

Также можно сделать так чтоб результат выполнения пайплайна не записывался в файл, а создавался NumPy array который вы можете использовать другими инструментами вашего воркфлоу или визуализировать в Jupyter.

Для этого нужно убрать из пайплайна этап записи результата в файл и добавить arrays = pipeline.arrays.

Ниже приведен минималистичный пример:

import json
import pdal

input_laz_file = "input.laz"

my_pipeline = [
    input_laz_file,
    {
        "type": "filters.expression",
        "expression": "Classification == 2"
    }
]

pipeline_json = json.dumps(my_pipeline)
pipeline = pdal.Pipeline(pipeline_json)
count = pipeline.execute()
arrays = pipeline.arrays
print(f"count: {count}")

Потом можно просто визуализировать в виде таблицы

import pandas as pd

arr = pipeline.arrays[0]
df = pd.DataFrame(arr)

# df.head()

# Если полей много, для предпросмотра можно выбрать только нужные:
df[["X", "Y", "Z", "Classification", "Red", "Green", "Blue", "Infrared"]].head()

Также можно просто визуализировать при помощи matplotlib

import numpy as np
import matplotlib.pyplot as plt

arr = pipeline.arrays[0]

# Для ускорения визуализации случайно выбираем индексы точек из массива arr
idx = np.random.choice(len(arr), size=min(50000, len(arr)), replace=False)

x = arr["X"][idx]
y = arr["Y"][idx]
z = arr["Z"][idx]

fig = plt.figure(figsize=(12, 7))
ax = fig.add_subplot(111, projection="3d")

sc = ax.scatter(
    x, y, z,
    c=z,          # красим по Z
    cmap="turbo", # палитра turbo
    s=0.2,
    vmin=z.min(),
    vmax=z.max()
)

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")

cbar = plt.colorbar(sc, ax=ax, pad=0.1)
cbar.set_label("Z")

plt.show()

Простые но полезные пайплайны

Ниже приведу несколько простых но полезных пайплайнов.

Обрезка облака точек по экстенту

Переменные minx, maxx, miny, maxy должны быть в той же системе координат что и исходный LAZ-файл.

pipeline_clip_by_extent = [
    {
        "type": "readers.las",
        "filename": input_laz_file
    },
    {
        "type": "filters.crop",
        "bounds": f"([{minx},{maxx}],[{miny},{maxy}])"
    },
    {
        "type": "writers.las",
        "filename": clipped_by_bbox_file
    }
]

Обрезка облака точек по полигону

PDAL принимает полигон в WKT формате. По этому нужно получить геометрию в WKT формате из исходного полигонального векторного файла. Для этого будем использовать geopandas.

import json
import pdal
import geopandas as gpd

def get_wkt_from_gpkg(gpkg_path):
    """Читает первый полигон из gpkg и возвращает его в формате WKT"""
    gdf = gpd.read_file(gpkg_path)
    # Берем геометрию первого объекта и объединяем, если там несколько полигонов
    geom = gdf.geometry.union_all()
    return geom.wkt

mask_wkt_geometry = get_wkt_from_gpkg(mask_file)

pipeline_clip_by_polygon = [
    input_laz_file,
    {
        "type": "filters.crop",
        "polygon": f"{mask_wkt_geometry}"
    },
    {
        "type": "writers.las",
        "filename": clipped_by_polygon_file
    }
]

pipeline_json = json.dumps(pipeline_clip_by_polygon)
pipeline = pdal.Pipeline(pipeline_json)
count = pipeline.execute()

Merging нескольких point clouds

Для объединения нескольких облаков точек нужно использовать PDAL инструмент filters.merge. Сначала в словарь передаем в pipeline список пути к нескольким LAZ файлам, и потом применяем фильтр filters.merge. На выходе получаем одно объединенное облако точек.

pipeline_merge = [
    input_laz_file_1,
    input_laz_file_2,
    {
        "type": "filters.merge"
    },
    {
        "type": "writers.las",
        "filename": output_merged_laz_file
    }
]

Создание GeoTIFF DEM из облака точек

DEM получается все дырявое. Дырки на месте вырезанных деревьев, зданий и других объектов. Дырки в растре можно потом заполнить другими инструментами.

target_resolution = 0.25  # Размер пикселя растра в единицах координат. В данном случае в метрах.

work_srs = "EPSG:3763"  # Система координат исходного облака точек.

pipeline_dem = [
    {
        "type": "readers.las",
        "filename": input_laz_2_dem_file
    },
    {
        "type": "filters.expression",
        "expression": "Classification == 2"
    },
    {
        "type": "writers.gdal",
        "filename": output_dem_geotiff_file,
        "override_srs": work_srs,  # необязательный параметр. Явно записывакет систему координат в GeoTIFF
        "resolution": target_resolution,
        "radius": target_resolution * 1.5,  # Радиус поиска точек вокруг центра пикселя.
        # радиус чуть больше размера пикселя помогает подобрать соседние точки даже там, где точки лежат не идеально по центрам пикселей.
        # слишком маленький radius → много пустых ячеек;
        # слишком большой radius → поверхность становится чрезмерно сглаженной.
        "output_type": "idw",  # способ, которым из соседних точек вычисляется значение пикселя.
        # "output_type" может быть: min, max, mean, idw, count, stdev
        "bounds": f"([{minx},{maxx}],[{miny},{maxy}])",  # границы выходного растра. Не обязательный параметр. Полезно чтобы точно контролировать extent итогового GeoTIFF
        "gdaldriver": "GTiff",
        "nodata": -9999  # Значение, которое будет записано в ячейки без данных
    }
]

Вы можете заполнять дырки со значением "no data" с помощью QGIS, GDAL, Rasterio, scipy и других инструментов.

Визуализация полученной цифровой модели рельефа в QGIS

Подведение итогов

Инструментов которые можно использовать в PDAL пайплайнах огромное количество. Пайплайны можно очень гибко настраивать и оптимизировать.

Надеюсь что после прочтения этой статьи вам будет проще и понятней разбираться с официальной документацией PDAL.

Устанавливать PDAL можно по инструкции с официального сайта PDAL quickstart.

Как я устанавливал PDAL, QGIS и некоторые другие спецефические библиотеки в Conda среду я описал в этом Telegram посте и комментариях к нему.