04 сентября 2011

Common Lisp. GIS системы.


PostGIS


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

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

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

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

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

PostGIS создает новый тип столбцов для Postgresql - geometry. Данный тип внутри себя может содержать следующие геометрические объекты:

  • POINT (точка)
  • LINESTRING (линия)
  • POLYGON (полигон)
  • MULTIPOINT (много точек)
  • MULTILINESTRING (много линий)
  • MULTIPOLYGON (много полигонов)
  • GEOMETRYCOLLECTION (много всего)
Сами данные вышеназванных типов внутри Postgresql хранятся внутреннем бинарном формате, и мы нуждаемся в некоторых методах вставки и получения геометрических объектов в и из базы. В стандарте OpenGIS для этого предусмотрены два вида функций: для текстового и бинарного (не такого как внутри Postgresql) представления объектов.
Текстовое представление служит для взаимодействия человека с машиной, бинарное, соответственно, участие человека не предполагает, а ориентируется на взаимодействие с другой программой. Названия данных представлений well-known text и well-known binary, (далее WKT и WKB).

PostGIS функции 


Ввод/вывод



Для ввода/вывода используются следующие SQL функции (они, повторюсь, заявлены в OpenGIS стандарте):
bytea WKB = ST_AsBinary(geometry);
text WKT = ST_AsText(geometry);

geometry = ST_GeomFromWKB(bytea WKB, SRID);
geometry = ST_GeometryFromText(text WKT, SRID);

Текстовое представление для типов объектов может выглядть, например, так:
  • POINT(0 0)
  • LINESTRING(0 0,1 1,1 2)
  • POLYGON( (0 0,4 0,4 4,0 4,0 0), (1 1, 2 1, 2 2, 1 2,1 1) )
  • MULTIPOINT(0 0,1 2)
  • MULTILINESTRING( (0 0,1 1,1 2), (2 3,3 2,5 4) )
  • MULTIPOLYGON( ( (0 0,4 0,4 4,0 4,0 0), (1 1,2 1,2 2,1 2,1 1) ), ( (-1 -1,-1 -2,-2 -2,-2 -1,-1 -1) ) )
  • GEOMETRYCOLLECTION(POINT(2 3), LINESTRING(2 3,3 4) )
С бинарным представлением чуть позже мы будем работать с помощью библиотеки cl-ewkb.

Кроме того PostGIS поддерживает текстовые вывод в таких форматах как:
  • GeoJSON
  • GML
  • KML
  • SVG
  • GeoHash

Сравнения


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

Операторы:

&& - Возвращает 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.

Функции-аналоги можно посмотреть здесь, но они в рамках статьи не понядобятся:

http://postgis.org/documentation/manual-1.5/reference.html#Spatial_Relationships_Measurements.



Преобразования



Функции можно посмотреть здесь: http://postgis.org/documentation/manual-1.5/reference.html#Geometry_Editors.

Из них нам могло бы понадобиться две функции для перемещения и для масштабирования, но кто-то из программистов PostGIS подумал за нас и предоставил одну большую функцию:


geometry ST_TransScale(geometry geomA, float deltaX, float deltaY, float XFactor, float YFactor);

Переносит объект на deltaX по X, на deltaY по Y, и умножает(масштабирует) на XFactor значения абсциссы, на YFactor значения ординаты.


Индексы



Индексы служат интструментом ускорения поиска. В нашем случае используются GiST тип индексов. Данный тип индексов используется только для функций/операторов сравнения, которые используют ограничивающие рамки. Например, индекс будет использоваться для оператора &&, и не будет при вычислении отношения: длина LINESTRING > 1000.

Данной информации нам пока достаточно.


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


Нужно было бы начать со скучной длительной истории о развертывании 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.

osm2pgsql

Эскпортирование из 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.

Рисование


Нам понадобится 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))

Вот они таблицы:
f_table_catalog f_table_schema f_table_name f_geometry_column coord_dimension srid type
public all_bounds the_geom 2 4326 MULTIPOLYGON
public osm_point way 2 4326 POINT
public osm_line way 2 4326 LINESTRING
public osm_polygon way 2 4326 GEOMETRY
public osm_roads way 2 4326 LINESTRING
Я хотел бы обратить внимание на проделанную Марижном работу. Пакет 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)
....

В анонимной функции нам необходимо создать новый контур, в месте его создания нарисуем подпись:
....
               (map 'vector (lambda (object) 
....
                    objects)
....

Функция 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 строчек кода:)

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

2 комментария:

  1. Классно! Ждём-с ещё статьи и переводы:)

    P.S. А ссылки на переводы не работают :(

    ОтветитьУдалить