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

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

== Введение ==
Лидар (LIDAR - Light Identification, Detection and Ranging) – это технология получения и обработки информации дистанционного зондирования с проведением высокоточных измерений X, Y, Z координат. Обработанные и пространственно организованные данные лидарной съемки называют «облаком точек» – это огромные наборы высотных 3D точек, имеющих дополнительную атрибутику.<br />
Данные LIDAR на территории городов или регионов могут содержать в себе миллиарды точек. Хранение их в базе данных в классическом виде как тип геометрии «точка» представляется возможным, но очень ресурсозатратным. Во-первых, огромное количество записей в таблице, что увеличивает время обработки данных. Во-вторых, большой объем памяти, занятый хранилищем.<br />
Для решения этих проблем существует технология хранения данных лидарной съемки в PostgreSQL/PostGIS c использованием встроенного расширения pointcloud и импорта таких данных в базу с помощью библиотеки PDAL (Point Data Abstract Library)<ref name="ttttest">[http://www.pdal.io www.pdal.io], Официальная страница PDAL.</ref>.<br />
В статье использованы материалы лидарной съемки на территорию Нижнего Новгорода, представленные файлами в формате LAS (Logging ASCII Standard) на площади размером 1×1 км.<br />
<br />
[[Файл:Lidar-data-nn-1.png|600px]]<br />
== Создание хранилища ==
Установка PostgreSQL/PostGIS для Windows описана на странице http://gis-lab.info/qa/postgis-install.html#02.<br />
Прежде чем загрузить данные LIDAR в таблицу PostgreSQL, а затем сделать пространственный анализ с помощью PostGIS, необходимо создать новую базу данных с соответствующими расширениями:<br />
# Создать новую базу данных с именем «Lidar».
# Активировать расширения pointcloud, postgis и pointcloud_postgis<ref name="tttest">[http://s3.cleverelephant.ca/foss4gna2013-pointcloud.pdf s3.cleverelephant.ca], LIDAR in PostgreSQL with PointCloud.</ref>.
<syntaxhighlight lang="sql">
CREATE EXTENSION postgis;
CREATE EXTENSION pointcloud;
CREATE EXTENSION pointcloud_postgis;</syntaxhighlight>
''Примечание: в ОС Windows расширения pointcloud и pointcloud_postgis устанавливается по умолчанию с PostGIS (http://postgis.net/windows_downloads/). Дополнительные настройки не требуются. При возникновении ошибок возможна отдельная загрузка расширения(https://github.com/pgpointcloud/pointcloud).''<br />
<br />
Объектами в pointcloud для PostGIS являются точки (PcPoint) - базовый объект, и блоки (PcPatch) - объединенные в группы PcPoint. Вместо таблицы миллиардов отдельных записей PcPoint, набор данных LIDAR может быть представлен в базе как гораздо меньшая коллекция (10 миллионов) из PcPatch записей.<br />
Расширение pointcloud_postgis добавляет функции, которые позволяют выполнять обработку облака точек, преобразовывать PcPoint и PcPatch в геометрию и осуществлять пространственную фильтрацию данных. Примеры использования таких функций рассмотрены в главе ''Пространственный анализ''. 
== Импорт данных ==
Импорт данных происходит с помощью библиотеки PDAL. PDAL это библиотека для манипулирования данными облака точек. Она очень похожа на библиотеку GDAL, которая обрабатывает растровые и векторные данные. В дополнение к коду библиотеки, PDAL предоставляет набор приложений командной строки, которые пользователи могут удобно использовать для обработки облака точек. <br />
=== Установка PDAL ===
Установить PDAL можно через OSGeo4W<ref name="ttttestt">[https://trac.osgeo.org/osgeo4w/ trac.osgeo.org/osgeo4w], Официальная страница OSGeo4W.</ref>.<br />
Для этого нужно скачать установщик с официального сайта OSGeo4W - https://trac.osgeo.org/osgeo4w/ и запустить его. Выбрать '''расширенную установку''' и следовать указаниям установщика. На вкладке «Выберите пакеты» в строке поиска ввести ''PDAL'' и начать установку.<br />
<br />
[[Файл:Window-install-pdal.png|600px]]<br />
<br />
Для проверки корректной работы PDAL, после установки, можно прочитать метаданные LAS-файла, выполнив команду:
<pre>pdal info --input H:\LIDAR\I_+1_+2_10.las --schema</pre>
Результат выполнения команды показан на рисунке. Команда –schema показывает размеры и типы полей: X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, ScanDirectionFlag, EdgeOfFlightLine, Classification, ScanAngleRank, UserData, PointSourceId и Time.<br />
<br />
[[Файл:Metadata-las-2_2.png]]<br />
<br />
=== Импорт данных в PostgreSQL через PDAL ===
Теперь необходимо создать файл загрузки для чтения LAS-файла, разбиения облака точек на блоки (например, по 400 точек), а затем записи блоков в базу данных.<br />
Файлом загрузки PDAL является XML-файл, который описывает обработку данных. Для загрузки нам понадобятся три драйвера PDAL:
*readers.las – для чтения LAS-файла
*filters.chipper - для разбиения облака точек на блоки
*writers.pgpointcloud – для записи LAS-файла в базу данных 
Создаем новый текстовый документ с именем «las2pg.xml» и расширением XML следующего вида<ref name="test">[http://workshops.boundlessgeo.com/tutorial-lidar/#loading-lidar-into-the-database workshops.boundlessgeo.com], Analyzing and Visualizing LIDAR.</ref>:
<syntaxhighlight lang="xml">
<?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></syntaxhighlight>
''Примечание: названия драйверов могут отличаться. Для своей версии PDAL узнать доступные драйвера можно с помощью команды pdal --drivers.''<br />
<br />
Теперь мы можем загрузить LAS-файл в базу данных. Для этого необходимо выполнить команду:
<pre>pdal pipeline H:\LIDAR\las2pg.xml</pre>
После завершения процесса загрузки в базе данных «LIDAR» появится новая таблица «Medford» следующей структуры:
<pre>                      Table "public.medford"
 Column |    Type    |                      Modifiers
--------+------------+------------------------------------------------------
 id     | integer    | default nextval('medford_id_seq'::regclass)
 pa     | pcpatch    |
 Indexes:
    "medford_pkey" PRIMARY KEY(id)</pre>
Хоть таблица и состоит из двух полей ID и PA, но в поле PA закодирована информация о точках лидара: X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, Classification, ScanAngleRank, Red, Green, Blue. В одной записи хранится пространственная и атрибутивная информация четырехсот точек. <br />
Для визуализации результата импорта можно воспользоваться любой ГИС с возможностью подключения к базе данных. Мы используем QGIS<ref name="teeest">[http://www.qgis.org/ru/site/ www.qgis.org], Официальная страница QGIS.</ref>, где создано соединение с базой данных «Lidar».<br />
<br />
[[Файл:Look-patch-lidar-qgis-4.png‎|600px]]<br />
<br />
Результатом такой технологии импорта данных LIDAR является существенное уменьшение размерности, по сравнению с обычным импортом графических слоев. Характеристики сравнения представлены в таблице.<br />
{| class="wikitable"
|-
! Показатель !! Импорт в 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.<br />
Для обработки и анализа импортированных данных мы можем использовать встроенные функции в расширении pointcloud<ref name="test" />. Далее приведены примеры использования таких функций.
<syntaxhighlight lang="sql">
-- Сколько точек в облаке (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</syntaxhighlight>
''Примечание: более подробную информацию о расширении pointcloud и функциях SQL смотри в документации (https://github.com/pgpointcloud/pointcloud, http://suite.opengeo.org/docs/latest/dataadmin/pointcloud/postgis.html).''
<br />
<br />
== Ссылки ==
<references />