Использование 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:

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

CREATE TEMPORARY TABLE w AS 
SELECT ST_Transfrom(l.geom, 28403) AS combine 
FROM water_line AS l
UNION ALL SELECT ST_Transfrom(ST_SetSRID(p.the_geom, 3395), 28403) 
FROM water_poly AS p;

CREATE INDEX temp_water_geom
  ON w
  USING gist
  (combine);

CREATE TEMPORARY TABLE a AS
SELECT DISTINCT ON (dump.gid) dump.gid,
  ST_Distance(ST_Transform(ST_SetSRID(dump.the_geom, 3395), 28403), w.combine) 
  AS buil_dist
FROM illegal_dumps AS dump  
LEFT JOIN w
FROM ST_Dwithin(ST_Transform(ST_SetSRID(dump.the_geom, 3395), 28403), w.combine), 10000)
WHERE dump."NAME"='illegal dump'
ORDER BY dump.gid, ST_Distance(ST_Transform(ST_SetSRID(dump.the_geom, 3395), 28403), w.combine);

UPDATE dump_dist
SET water_dist = (SELECT CAST(buil_dist AS int) FROM a WHERE dump_dist.gid = a.gid)
FROM a
WHERE dump_dist.gid = a.gid;

DROP TABLE a, w;