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

Введение[править]

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

Создание хранилища[править]

Установка PostgreSQL/PostGIS для Windows описана на странице http://gis-lab.info/qa/postgis-install.html#02.
Прежде чем загрузить данные LIDAR в таблицу PostgreSQL, а затем сделать пространственный анализ с помощью PostGIS, необходимо создать новую базу данных с соответствующими расширениями:

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

CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;

Примечание: в ОС Windows расширения pointcloud и pointcloud_postgis устанавливается по умолчанию с PostGIS (http://postgis.net/windows_downloads/). Дополнительные настройки не требуются. При возникновении ошибок возможна отдельная загрузка расширения(https://github.com/pgpointcloud/pointcloud).

Объектами в pointcloud для PostGIS являются точки (PcPoint) - базовый объект, и блоки (PcPatch) - объединенные в группы PcPoint. Вместо таблицы миллиардов отдельных записей PcPoint, набор данных LIDAR может быть представлен в базе как гораздо меньшая коллекция (10 миллионов) из PcPatch записей.
Расширение pointcloud_postgis добавляет функции, которые позволяют выполнять обработку облака точек, преобразовывать PcPoint и PcPatch в геометрию и осуществлять пространственную фильтрацию данных. Примеры использования таких функций рассмотрены в главе Пространственный анализ.

Импорт данных[править]

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

Установка PDAL[править]

Установить PDAL можно через OSGeo4W^[3].
Для этого нужно скачать установщик с официального сайта OSGeo4W - https://trac.osgeo.org/osgeo4w/ и запустить его. Выбрать расширенную установку и следовать указаниям установщика. На вкладке «Выберите пакеты» в строке поиска ввести PDAL и начать установку.



Для проверки корректной работы 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.

Импорт данных в PostgreSQL через PDAL[править]

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

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

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

<?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^[5], где создано соединение с базой данных «Lidar».



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

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

Пространственный анализ[править]

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

-- Сколько точек в облаке (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

Примечание: более подробную информацию о расширении pointcloud и функциях SQL смотри в документации (https://github.com/pgpointcloud/pointcloud, http://suite.opengeo.org/docs/latest/dataadmin/pointcloud/postgis.html).

Ссылки[править]