Импорт и анализ данных лидарной съемки в PostgreSQL/PostGIS с помощью PDAL

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Эта страница является черновиком статьи.


Рассмотрен процесс импорта и анализа данных лидарной съемки в PostgreSQL с использованием технологии Point Cloud и библиотеки PDAL

Введение

Лидар (LIDAR - Light Identification, Detection and Ranging) – это технология получения и обработки информации дистанционного зондирования с проведением высокоточных измерений X, Y, Z координат. Обработанные и пространственно организованные данные лидарной съемки называют облаком точек – это огромные наборы высотных 3D точек, имеющих дополнительную атрибутику.
Данные LIDAR на территории городов или регионов могут содержать в себе миллиарды точек местности. Хранение их в базе данных в классическом виде как тип геометрии «точка» представляется возможным, но очень ресурсозатратным. Во-первых, огромное количество записей в таблице, что увеличивает время обработки данных. Во-вторых, большой объем памяти, занятый хранилищем.
Для решения этих проблем существует технология хранения и обработки данных лидарной съемки в PostgreSQL/PostGIS c помощью библиотеки PDAL (Point Data Abstract Library).
В статье использованы материалы лидарной съемки на территорию Нижнего Новгорода, представленные файлами в формате LAS (Logging ASCII Standard) на площади размером 1×1 км.

Lidar-data-nn-1.png

Создание хранилища

Прежде чем загрузить данные LIDAR в таблицу PostgreSQL, а затем сделать пространственный анализ с помощью PostGIS, необходимо создать новую базу данных с соответствующими расширениями:

  1. Создать новую базу данных с именем «Lidar».
  2. Создать расширения PointCloud, PostGIS и pointcloud_postgis[1].
CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;

Импорт данных

Импорт данных происходит с помощью библиотеки PDAL. PDAL это библиотека для манипулирования данными облака точек. Она очень похожа на библиотеку GDAL, которая обрабатывает растровые и векторные данные. В дополнение к коду библиотеки, PDAL предоставляет набор приложений командной строки, которые пользователи могут удобно использовать для обработки облака точек. Установить PDAL можно через OSGeo4W.
Для проверки корректной работы PDAL, после установки, можно прочитать метаданные LAS-файла, выполнив команду: 

pdal info --input H:\LIDAR\I_+1_+2_10.las --schema

Результат выполнения команды показан на рисунке. Команда –schema показывает размеры и типы полей: X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, ScanDirectionFlag, EdgeOfFlightLine, Classification, ScanAngleRank, UserData, PointSourceId и Time.

Metadata-las-2 2.png

Теперь необходимо создать файл загрузки для чтения LAS-файла, разбиения облака точек на блоки (например, по 400 точек), а затем записи блоков в базу данных.
Файлом загрузки PDAL является XML-файл, который описывает обработку данных. Для загрузки нам понадобятся три драйвера PDAL:

  • writers.pgpointcloud – для чтения LAS-файла
  • filters.chipper
  • readers.las – для записи LAS-файла в базу данных

Создаем новый текстовый документ с именем «las2pg.xml» и расширением XML следующего вида[2]: 

<?xml version="1.0" encoding="utf-8"?>
<Pipeline version="1.0">
  <Writer type=" writers.pgpointcloud ">
    <Option name="connection">dbname=lidar user=postgres</Option>
    <Option name="table">medford</Option>
    <Option name="srid">4326</Option>
    <Filter type=" filters.chipper">
      <Option name="capacity">400</Option>
        <Reader type="readers.las">
          <Option name="filename"> H:\LIDAR\I_+1_+2_10.las</Option>
          <Option name="spatialreference">EPSG:4326</Option>
        </Reader>
    </Filter>
  </Writer>
</Pipeline>

Примечание: названия драйверов могут отличаться. Для своей версии PDAL узнать доступные драйвера можно с помощью команды pdal --drivers.
Теперь мы можем загрузить LAS-файл в базу данных. Для этого необходимо выполнить команду: 

pdal pipeline H:\LIDAR\las2pg.xml

После завершения процесса загрузки в базе данных «LIDAR» появится новая таблица «Medford» следующей структуры: 

                      Table "public.medford"
 Column |    Type    |                      Modifiers
--------+------------+------------------------------------------------------
 id     | integer    | default nextval('medford_id_seq'::regclass)
 pa     | pcpatch    |
 Indexes:
    "medford_pkey" PRIMARY KEY(id)

Хоть таблица и состоит из двух полей ID и PA, но в поле PA закодирована информация о точках лидара: X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, Classification, ScanAngleRank, Red, Green, Blue. В одной записи хранится пространственная и атрибутивная информация четырехсот точек.
Для визуализации результата импорта можно воспользоваться любой ГИС с возможностью подключения к базе данных. Мы используем QGIS, где создано соединение с базой данных «Lidar».

Look-patch-lidar-qgis-4.png

Результатом такой технологии импорта данных LIDAR является существенное уменьшение размерности, по сравнению с обычным импортом графических слоев. Характеристики сравнения представлены в таблице.

Показатель Импорт в PostgreSQL/PostGIS через PDAL Импорт в PostgreSQL/PostGIS LAS-файла
Количество записей 1 864 931 919
Размер таблицы, MB 9,9 197

Таким образом импортирован один LAS-файл. Если таких файлов множество, есть два пути решения. Первый – объединить все LAS-файлы в один, например, с помощью набора программного обеспечения LAStools (LAStools(с) by Martin Isenburg), в частности программы las2las. Второй – написать собственный скрипт, который генерирует файлы загрузки и выполняет команду PDAL.

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

Формат вывода получившейся таблицы в PostgreSQL труден для восприятия, поэтому для пространственного анализа можно отформатировать таблицу с данными LIDAR.
Для обработки и анализа импортированных данных мы можем использовать встроенные функции в расширении PointCloud[2]. 

-- Сколько точек в облаке (931919)
SELECT Sum(PC_NumPoints(pa))
FROM medford;

-- Какова средняя высота первого блока? (151.95)
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford LIMIT 1
)
SELECT Avg(PC_Get(pt,'Z')) FROM pts;

-- Какой вид имеет первая точка?
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford LIMIT 1
)
SELECT PC_AsText(pt) FROM pts LIMIT 1;
-- {
--  "pcid":1,
--  "pt":[255,1,1,0,0,2,0,0,10,556976,2250.04,1253.42,151.95]
-- }

-- Сколько блоков? (1864)
SELECT Count(*)
FROM medford;

-- Какая минимальная и максимальная высота в блоке? (141.78/ 204.33)
SELECT
  Min(PC_PatchMin(pa, 'z')) AS min,
  Max(PC_PatchMax(pa, 'z')) AS max
FROM medford;
-- Как выглядит геометрия блока?
SELECT st_asewkt(pa::geometry) FROM medford LIMIT 1;
-- SRID=4326;POLYGON((2250.04 1250.2,2250.04 1253.82,2258.88 1253.82,
-- 2258.88 1250.2,2250.04 1250.2))

-- Как выглядит геометрия точки?
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford LIMIT 1
)
SELECT ST_AsEWKT(pt::geometry) FROM pts LIMIT 1;
-- SRID=4326;POINT(2250.04 1253.42 151.95)

-- Выбрать атрибуты точек: тип классификации, координаты Z,X,Y
WITH pts AS (
  SELECT PC_Explode(pa) AS pt
  FROM medford 
)
SELECT PC_Get(pt,'Classification'),PC_Get(pt,'Z'), PC_Get(pt,'X'), PC_Get(pt,'Y') FROM pts LIMIT 1;
-- 2
-- 151.95
-- 2250.04
-- 1253.42

Примечание: более подробную информацию о расширении Point Cloud и функциях SQL смотри в документации[3].

Ссылки по теме:

  1. ↑ s3.cleverelephant.ca, LIDAR in PostgreSQL with PointCloud.
  2. ↑ 2,0 2,1 workshops.boundlessgeo.com, Analyzing and Visualizing LIDAR.
  3. ↑ github.com, A PostgreSQL extension for storing point cloud (LIDAR) data.
Источник — https://wiki.gis-lab.info/index.php?title=Импорт_и_анализ_данных_лидарной_съемки_в_PostgreSQL/PostGIS_с_помощью_PDAL&oldid=24588