Использование PostGIS для измерения расстояний

Эта страница является черновиком статьи.

При работе с пространственными данными часто возникает необходимость вычислить расстояния от одних объектов до других. Для решения этой задачи как нельзя лучше подходит PostGIS.

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

Задача

Вычислить расстояния от некоторых объектов данного слоя до ближайших водоёмов или рек.

Данные

В нашей базе данных содержатся следующие таблицы:

Название таблицы	Тип геометрии	Проекция	Первичный ключ	Примечания
water_line	полилинии	EPSG: 3395	gid	очень много объектов
water_poly	мультиполигоны	не задана (EPSG: 3395)	gid	очень много объектов
illegal_dump	мультиполигоны	не задана (EPSG: 3395)	gid	мало объектов
dump_dist	без поля геометрии	без проекции	gid	предназначена для хранения расстояний от свалок до различных объектов

Как мы видим, все используемые таблицы, содержащие поле геометрии, находятся в одной проекции (точнее это мы знаем, что они в одной проекции), однако, в самой базе данных, проекция для water_poly и illegal_dump не задана. Кроме того, используемая проекция не подходит для измерения расстояний. Таблицы с водными объектами содержат десятки и даже сотни тысяч записей, что делает вычисление расстояний до них весьма продолжительным. Ну и напоследок самое прекрасное: нас интересует расстояние до любого ближайшего водного объекта и неважно, находится ли он в таблице water_line или water_poly.

Планирование запроса

Очевидно, что при расчёте расстояний, нам понадобиться задать проекцию таблицам, где она не задана; осуществить перепроецирование в проекцию более подходящую для расчёта расстояний (соответствующую зону Гаусса-Крюгера); совместить геометрии полилиний и мультиполигонов в одном поле геометрии (да, в отличие от .shp, PostGIS позволяет хранить различные типы геометрий в одном поле); ну и наконец, заполонить таблицу расстояний искомыми значениями.

Разобьём планируемый запрос на несколько логических этапов:

1. Объединение и перепроецирование геометрий водных объектов из разных таблиц.
2. Создание пространственного индекса для объединённых и перепроецированных геометрий водных объектов (для последующего ускорения запроса для вычисления расстояний).
3. Непосредственно вычисление расстояний от определённых объектов таблицы illegal_dumps до ближайшего водного объекта.
4. Запись полученных значений в таблицу расстояний (уточним, что в таблице dump_dist поля имеют тип данных int, чтобы хранить расстояния в целых метрах в то время, как функция определения расстояния возвращает тип float8).
5. Удаление временных файлов.

Нам понадобятся следующие функции PostGIS/PostgreSQL:

• ST_Distance
• ST_Transfrom
• ST_SetSRID
• ST_DWithin
• SELECT
• CREATE TABLE
• DROP TABLE
• CREATE INDEX
• UNION ALL
• LEFT JOIN
• ORDER BY
• UPDATE
• CAST

Составление запроса

Для начала создадим временную таблицу, в которую поместим перепроецированные геометрии из таблиц, содержащих водные объекты:

CREATE TEMPORARY TABLE w AS  -- создаём временную таблицу
SELECT ST_Transfrom(l.geom, 28403) -- выбираем все записи поля геометрии таблицы water_line и перепроецируем их в EPSG:28403
  AS combine -- задаём имя результирующего поля геометрии 
FROM water_line AS l -- задаём псевдоним таблицы для данного запроса
UNION ALL -- объединяем полученный результат со следующим запросом:
  SELECT ST_Transfrom(ST_SetSRID(p.geom, 3395), 28403) -- выбираем все записи поля геометрии из water_poly,
                                                       -- задаём им проекцию EPSG:3395 и перепроецируем в EPSG:28403
FROM water_poly AS p; -- задаём псевдоним таблицы для данного запроса

Для того, чтобы существенно ускорить последующий запрос для вычисления расстояний, создадим пространственный индекс для нашей временной таблицы w:

CREATE INDEX temp_water_geom -- создаём индекс и даём ему имя
  ON w -- определяем индексируемую таблицу
  USING gist -- указываем тип индекса
  (combine); -- указываем индексируемое поле

Опционально можно создать перманентный пространственный индекс для геометрии таблицы illegal_dump в проекции 28403:

CREATE INDEX pulkovo_dumps
  ON illegal_dumps
  USING gist
  (ST_Transform(ST_SetSRID(dump.geom, 3395), 28403));

Теперь вычислим расстояния. Нас интересуют те объекты таблицы dump_dist, которые в поле NAME имеют значение "illegal dump". Чтобы существенно ускорить процесс (PostGIS будет вычислять расстояния от каждого объекта до каждого объекта, что может занять очень много времени) отсеем те водные объекты, которые, находятся на расстоянии, скажем 10000 метров от свалок (мы уверены, что расстояние от свалок до водных объектов не превышает 10 км). Для этого применим функцию ST_DWithin, которая и воспользуется пространственным индексом, который мы создали для временной таблицы w.

CREATE TEMPORARY TABLE a AS
SELECT DISTINCT ON (dump.gid) dump.gid, -- задаём условие, чтобы одной записи в таблице illegal_dumps 
                                        -- соответствовало только одно значение расстояния
  ST_Distance(ST_Transform(ST_SetSRID(dump.geom, 3395), 28403), w.combine) -- вычисляем расстояния между объектами двух таблиц
  AS wat_dist
  FROM illegal_dumps AS dump  
LEFT JOIN w
  FROM ST_Dwithin(ST_Transform(ST_SetSRID(dump.geom, 3395), 28403), w.combine), 10000) -- убираем из расчётов 
                           -- все объекты таблицы w, которые находятся дальше 10 км от объектов illegal_dumps
  WHERE dump."NAME"='illegal dump'
ORDER BY dump.gid, -- сортируем результата сначала по значениям gid, 
         ST_Distance(ST_Transform(ST_SetSRID(dump.geom, 3395), 28403), w.combine); --а затем - по значениям расстояний 
                                                                                   -- между объектами (таким образом 
                                                 -- в результирующей таблице остаются только самые короткие расстояния)

Теперь занесём полученные значения в таблицу dump_dist (поле water_dist). Для того, чтобы преобразовать значения типа float8 (нам ведь не нужна нанометровая точность), которые нам вернула функция ST_Distance, в int (тип данных поля water_dist - int) используем функцию CAST.

UPDATE dump_dist -- обновляем таблицу dump_dist
SET water_dist =  -- устанавливаем значения поля water_dist равными:
               (SELECT CAST(buil_dist AS int) -- выбираем значения из таблицы "a" и задаём им тип int.
                FROM a) 
FROM a
WHERE dump_dist.gid = a.gid;

Удаляем ставшие более ненужными временные таблицы a и w.

DROP TABLE a, w;

Заключение

На простом примере из реальной практики было продемонстрировано использование комбинирование различных функций PostgreSQL/PostGIS для расчёта расстояний между объектами.