PostGIS

В продолжение темы о Postgresql. Я думаю многим известно, что в этой чудесной СУБД, помимо бумажек, кадров, прибылей и убылей, можно хранить геометрические данные. Для этого достаточно установить расширение PostGIS.

Данное расширение предоставляет:

• геометрические типы данных
• функции для работы с ними
• таблицы с метаданными
• индексы

Кроме того, PostGIS реализует интерфейс OpenGIS, что позволяет создавать приложения, не зависящие от реализации гео базы данных.

Рассмотрим его слегка.

bytea WKB = ST_AsBinary(geometry);

text WKT = ST_AsText(geometry);

geometry = ST_GeomFromWKB(bytea WKB, SRID);

geometry = ST_GeometryFromText(text WKT, SRID);

&& - Возвращает TRUE если ограничивающая рамка A пересекается c рамкой B.
&< - Возвращает TRUE если ограничивающая рамка A пересекается или левее рамки B.
&<| - Возвращает TRUE если ограничивающая рамка A пересекается или ниже рамки B.
&> - Возвращает TRUE если ограничивающая рамка A пересекается или правее рамки B.
<< - Возвращает TRUE если ограничивающая рамка A строго слева от рамки B.
<<| - Возвращает TRUE если ограничивающая рамка A строго ниже рамки B.
= - Возвращает TRUE если ограничивающая рамка A совпадает с рамкой B.
>> - Возвращает TRUE если ограничивающая рамка A строго правее рамки B.
@ - Возвращает TRUE если ограничивающая рамка A содержится в рамке B.
|&> - Возвращает TRUE если A's ограничивающая рамка пересекается или выше рамки B.
|>> - Возвращает TRUE если A's ограничивающая рамка строго выше рамки B.
~ - Возвращает TRUE если ограничивающая рамка A содержит рамку B.
~= - Возвращает TRUE если A's ограничивающая рамка совпадает с рамкой B.

Развертывание

Нужно было бы начать со скучной длительной истории о развертывании PostGIS'а, но мне повезло. Я натолкнулся на гостевой доступ к базе данных, которая, кроме вышеупомянутого расширения, еще и содержит географические данные стран бСССР. Данные предоставлены проектом OpenStreetMap (далее OSM). OSM некоторое время назад мигрировал с MySql на Postgres и, еще не успев применить PostGIS, хранит данные как есть (колонки latitude и longitude). Gis-lab'овцы импортируют часть данных из OSM с помощью программы osm2pgsql, которая и создает PostGIS объекты.

Сначала расскажу про модель данных, используемую в OSM. Она очень простая.

OpenStreetMap

Ноды: Точки используемые для пометок определенных мест или для соединения сегментов.

Пути: Упорядоченный список нод для отображения сегментов линии. Используется для дорог, путей и т.д.

Закрытые пути (полигоны): закрытые пути - это закольцованная линия. Используется для отображения парков, озер, островов, зданий и т.д.

Отношения: Когда различные пути соединены между собой, но не представляют один и тот же физический объект, используется отношение для описании функции каждого из пути. Они используются для описания таких вещей, как велодорожки, "turn restrictions" и участки с отверстиями.

Данные объекты имеют теги. Теги выполняют описательную роль, используются при визуализации карты для стилизации объектов и создания подписей, а также при поиске объектов в базе данных.

Подробнее: http://wiki.openstreetmap.org/wiki/Beginners_Guide_1.3.

OpenStreetMap 2 PostGIS

Эскпортирование из OSM осуществляется в xml файл. Импортирование из данного файла в Postgresql/PostGIS базу данных выполняется программой osm2pgsql. Вот небольшое описание схемы импортирования. Взято здесь: http://wiki.openstreetmap.org/wiki/Osm2pgsql/schema.

Нижеперечисленные таблицы содержат географические данные:

• osm_line: содержит все импортированные пути.
• osm_polygon: содержит все импортированные полигоны.
• osm_point: содержит все импортированные ноды с тэгами.
• osm_roads: содержит подмножество 'osm_line' предназначенное для низкого разрешения. Выборка производится в соответствии с тэгами (какими не известно).

Каждая из таблиц имеет столбец way, который содержит геометрические данные. По два индекса созданы для каждой из таблиц: один на столбец osm_id и один на way. Координаты геометрических объектов в проекции EPSG:900913 AKA G00GlE. Примечание. На самом деле используется проекция указанная в столбце srid таблицы geometry_columns. В нашем случае это EPSG:4326. Просмотреть все возможные проекции можно в таблице spatial_ref_sys.

Отношения напрямую не импортируются, а представляют собой несколько строк в таблице osm_line.

gis-lab.info

Нам понадобится quicklisp, postmodern, cl-ewkb, и vecto. quicklisp скачайте и установите. Уже в slime разрешите зависимости.

(mapcar #'ql:quickload '(:postmodern :cl-ewkb :vecto))

Сервер баз данных запущен на gis-lab.info. Подключимся к нему:

(in-package :postmodern) 
(connect-toplevel "osm" "guest" "guest" "gis-lab.info")

Можете просмотреть список таблиц:

(list-tables)

Их будет много, не обращайте внимания. Нам понадобяться лишь несколько из них. Метаданные о "геометрических" столбцах хранятся в таблице geometry_columns:

(query (:select '* :from 'geometry_columns))

Так неудобно, не видно названий столбцов. Давайте так:

(query (:select '* :from 'geometry_columns) :plists)

Предыдущий вариант также не удобен, давайте так.

(defvar headers (mapcar (lambda (x) (car x)) (table-description 'geometry_columns)))
(defvar rows (query (:select '* :from 'geometry_columns)))
(format nil "~{ ~{ ~19< ~A ~> ~^|~} ~% ~}" (cons headers rows))

Вот они таблицы:

Я хотел бы обратить внимание на проделанную Марижном работу. Пакет S-SQL позволяет записать SQL запрос в терминах и синтаксисе лиспа. Keyword'ы представляют собой ключевые слова SQL, символы транслируются в SQL идентификаторы, ключевые слова в начале списков - в вызовы функций или применения операторов. Подробности вы можете узнать в переводе руководства к данному пакету: ./libraries%3Apostmodern#%D0%A1%D0%BF%D1%80%D0%B0%D0%B2%D0%BE%D1%87%D0%BD%D0%BE%D0%B5%20%D1%80%D1%83%D0%BA%D0%BE%D0%B2%D0%BE%D0%B4%D1%81%D1%82%D0%B2%D0%BE%20S-SQL. Одним из плюсов S-SQL является тот факт, s-выражения получаются структурированными, в отличии от строкового SQL запроса. Это позволяет автоматически форматировать их, и облегчает визуальное "проигрывание" кода.

Теперь давайте глянем на таблицу путей для большого масштаба поближе:

(format nil "~:{ ~A ~% ~}" (table-description 'osm_roads))
" osm_id 

  note

.....

  wood 

  way 

 "

Под многоточием понимается список столбцов. Для каждого тега для объекта, представленного в таблице, заведен отдельный столбец. Логическое значение тега можно просмотреть по ссылке http://wiki.openstreetmap.org/wiki/Key:TAGNAME. Например, граница государства имеет тег boundary со значением равным *adminstrative*. При этом граница именно государства иммет административный уровень (тег admin_level) равный двум. А границы подчиненных государству административных територий могут иметь admin_level от 3 до 10.

Как вы помните, индексы созданы только для столбцов *osm_id* и *way*. Засчет последнего индекса мы с легкостью может попросить данные, пересекающиеся с некоторым квадратом. Например, данные между 27 и 28 долготами и 54 и 55 широтами можно получить таким запросом.

;; SELECT way FROM osm_roads WHERE way && ST_GeometryFromText 'LINESTRING(27.0 54.0, 28.0 55.0)'
(query (:select 'way :from 'osm_roads :where (:&& 'way (:ST_GeometryFromText "LINESTRING(27.0 54.0, 28.0 55.0)"))))

Оператор пересечения && находит ограничивающие рамки левого и правого объектов и возвращает T, если рамки пересекаются. Мы могли бы использовать функцию ST_Intersects, но я хочу показать, как Postmodern S-SQL можно научить ранее неизвестным операторам. Сейчас он, понятное дело, не работает, так как модуль S-SQL ничего о нем не знает. Решение заключается в регистрации нового оператора:

(register-sql-operators :2+-ary :&&)

Вторым аргументом мы указали "арность" оператора, в данном случае оператор имеет смысл только при двух аргументах. Так как postmodern не имеет символа строгой 2-арности, используем "2 и более арность".

Мы использовали LINESTRING, который представляет диагональ квадрата, которым мы и охватываем необходимые данные.

Но не торопитесь, то, что мы получили - это просто внутреннее представление геометрического типа PostGIS. OpenGIS стандарт предполагает, что если мы хотим получить текстовый формат, мы должны преобразовать данные функцией ST_AsText. Но опять не торопитесь, у нас нет парсера текстового представления, есть только парсер бинарного в пакете cl-ewkb. Поэтому мы должны преобразовать результат SQL функцией ST_AsBinary, а затем уже в lisp-е полученный список отобразить распарсив well-known binary данные.

Теперь вопрос касаемый рисования. Данные мы получим в квадрате (27.0 54.0, 28.0 55.0). Но библиотека vecto позволяет рисовать объекты, координаты которых целые числа. Для этого мы можем на SQL стороне масштабировать изображение. Я предлагаю увеличить в 1000 раз, соответственно в будущем размер холста у нас будет также 1000x1000, а также предлагаю сместить полученные объекты на начало координат. Все это производится функцией ST_TransScale:

(defvar objects
            (mapcar (lambda (x)
                      (list (car x)
                            (cl-ewkb:decode (cadr x))))
                    (query (:select 'name 
                                               (:ST_AsBinary 
                                                (:ST_TransScale 'way -27.0 -54.0 1000 1000))
                                               :from 'osm_roads
                                               :where (:&& 'way
                                                           (:ST_GeometryFromText "LINESTRING(27.0 54.0, 28.0 55.0)"))))))

SQL уровень

• выбираем объекты, геометрия которых пересекаеться с квадратом bottom-left = 27.0 54.0 top-right = 28.0 55.0
• перемещаем вектором (-27.0 -54.0)
• увеличиваем вектором (1000 1000)
• отображаем объекты в well-known binary.

LISP уровень

• отображаем объекты из well-known binary в lisp струкутры.

Мы получаем столбец бинарных данных, пропускаем его через функцию декодирования cl-ewkb:decode и получаем вектор структур, которые нам надо отрисовать.

Рисование

Итак, у нас есть массив структур cl-ewkb:line-string, каждая из которых в свою очередь содержит массив структур cl-ewkb:point-primitive. Последняя содержит в себе координаты.

Сначала я хотел рисовать с помощью библиотеки cl-svg и даже начал переводить руководство. Данная библиотека позволяет создавать векторную графику в формате SVG. Данный формат - это обычный xml, и библиотека по сути генерирует текстовый файл. Но когда я столкнулся с необходимостью переноса/поворота системы координат для географического представления, энтузиазм резко пропал.

Тогда я обратился к библиотеке vecto. Что интересно, это pure-lisp решение включая зависимости, и у нас будет шанс оценить прикладную скорость реализации common lisp'а. Ну ее я собственно быстренько перевел, качество средненькое получилось. Данная библиотека позволяет рисовать из нижнего левого угла, умеет рисовать текст, ну и этого нам пока хватит.

Повторюсь, с символом objects у нас связан вектор геометирческих объектов и их подписей. Попробуем их отрисовать.

Экспортируем символ objects из postmodern, перейдем в пакет vecto, определим холст с помощью макроса with-canvas:

(export 'objects)
(in-package :vecto)
(defvar objects postmodern:objects)
(with-canvas (:width 1000 :height 1000)
....
)

Загрузим шрифт:

....

             (let ((font (get-font "/usr/share/fonts/TTF/times.ttf")))
....

Установим некоторые параметры рисования (цвет, ширину линии, размер шрифта):

....

               (set-rgb-stroke 1 0 0)

               (set-line-width 1)

               (set-font font 14)

....

Произведем итерацию по объектам полученным из базы данных:

....

               (map 'vector (lambda (object) 
....
                    objects)
....

В анонимной функции нам необходимо создать новый контур, в месте его создания нарисуем подпись:

....

                              (let ((point (elt (cl-ewkb:line-string-points-primitive (cadr object)) 0)))
                                (move-to 
                                 (round (cl-ewkb:point-primitive-x point))
                                 (round (cl-ewkb:point-primitive-y point)))
                                (when (postmodern:coalesce (car object))
                                  (draw-string 
                                   (round (cl-ewkb:point-primitive-x point))
                                   (round (cl-ewkb:point-primitive-y point)) (car object)))
                                )

.....

Функция postmodern:coalesce возвращает первый не-:NULL (именно keyword) аргумент. Если аргумент не найден она возвращает NULL.

В этой же анонимной функции рисуем все линию. Напомню, что в SQL запросе мы переместили и смасштабировали координты объектов, сейчас нам необходимо просто их округлить:

....

                              (map 'vector (lambda (point) 
                                             (line-to 
                                              (round (cl-ewkb:point-primitive-x point))
                                              (round (cl-ewkb:point-primitive-y point))))
                                   (cl-ewkb:line-string-points-primitive (cadr object)))

....

В ней же фиксируем результаты:

....

                              (stroke))

....

Теперь сохраняем результат в png файл:

....

               (save-png "test.png")))

Весь блок кода ответсвенный за отрисовку:

(export 'objects)
(in-package :vecto)
(defvar objects postmodern:objects)
(with-canvas (:width 1000 :height 1000)
             (let ((font (get-font "/usr/share/fonts/TTF/times.ttf")))
               (set-rgb-stroke 1 0 0)
               (set-line-width 1)
               (set-font font 16)
               (map 'vector (lambda (object) 
                              (let ((point (elt (cl-ewkb:line-string-points-primitive (cadr object)) 0)))
                                (move-to 
                                 (round (cl-ewkb:point-primitive-x point))
                                 (round (cl-ewkb:point-primitive-y point)))
                                (when (postmodern:coalesce (car object))
                                  (draw-string 
                                   (round (cl-ewkb:point-primitive-x point))
                                   (round (cl-ewkb:point-primitive-y point)) (car object)))
                                )
                              (map 'vector (lambda (point) 
                                             (line-to 
                                              (round (cl-ewkb:point-primitive-x point))
                                              (round (cl-ewkb:point-primitive-y point))))
                                   (cl-ewkb:line-string-points-primitive (cadr object)))
                              (stroke))
                    objects)
               (save-png "test.png")))

Заключение

Отмечу нереализованные задачи:

• Запрос только определенных гео объектов для разного масштаба.
• Запрос атрибутов объектов для стиля отрисовки.
• Атрибуты также влияют на z-order.
• Отрисовка масштабной линейки/сетки.
• Стилизованная отрисовка подписей к объектам.

То что получилось, конечно, представляет собой простой GIS helloworld, но несмотря на такой большой список задач, мне кажется, то, что есть - довольно неплохой результат: возможность нарисовать карту любой точки земли за 50 строчек кода:)

Пожелания, критика приветствуются.