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

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


Рассмотрен процесс хранения и анализа данных лидарной съемки в 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 км.

Lidar-data-nn-1.png

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

Установка 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 и начать установку.

Window-install-pdal.png

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

Импорт данных в 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».

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[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).

Ссылки

  1. ↑ www.pdal.io, Официальная страница PDAL.
  2. ↑ s3.cleverelephant.ca, LIDAR in PostgreSQL with PointCloud.
  3. ↑ trac.osgeo.org/osgeo4w, Официальная страница OSGeo4W.
  4. ↑ 4,0 4,1 workshops.boundlessgeo.com, Analyzing and Visualizing LIDAR.
  5. ↑ www.qgis.org, Официальная страница QGIS.
Источник — https://wiki.gis-lab.info/index.php?title=Импорт_и_анализ_данных_лидарной_съемки_в_PostgreSQL/PostGIS_с_помощью_PDAL&oldid=24895