Рисуем на тайлах электронной карты в MSSQL
Хочу рассказать читателям хабра-сообщества как используя CLR библиотеку Microsoft.SqlServer.Types можно формировать тайлы для электронной карты. В статье пойдёт речь о генерации списка картографических тайлов для их дальнейшего рендеринга. Будет описан алгоритм генерации тайлов по геометрии объектов, хранящейся в базе данных MS SQL 2008. Весь процесс рендеринга шаг за шагом будет рассмотрен на примере в конце статьи.
Содержание
Проблема
Исходные данные
Решение
Хранилище тайлов
Этапы подготовки тайлов
Используемые функции
Пример с ломаной линией
Проверка пересечения
Таблицы для хранения образов тайлов
Размещение иконки на тайле
Объединение тайлов
Отрисовка геометрии на тайле
Заключение
Проблема
Когда в браузере отображается большое количество гео данных в векторной графике (с помощью SVG или CANVAS), иногда приходится ждать не только пока данные загрузятся на клиентскую машину, но и пока выполнится процесс отрисовки, который может выполняться слишком долго.
При отображении большого количества иконок на карте в браузере можно применить кластеризацию, но для сложных геометрических объектов нужно использовать другой подход.
Исходные данные:
Набор геометрических объектов хранится в таблице базы данных Microsoft SQL 2008. Координаты узлов – это широта и долгота (EPSG:4326). Поле с гео-данными имеет тип GEOMETRY. Объекты должны отображаться на карте в виде иконки, для геометрии Point. В виде ломаной линии определённой толщины для геометрии Polyline. Геометрия Polygon должны отображаться в виде одного или нескольких закрашенных многоугольников с контуром. Тайлы должны соответствовать проекции Web Mercator
Решение:
Вместо векторной графики будем отображать объекты на карте в виде растрового слоя, то есть такими же изображениями (тайлами) как и сама карта. Для этого необходимо подготовить набор картографических тайлов для каждого масштаба с изображением объектов. Для формирования тайлов будем использовать проекцию Google Web Mercator, то есть преобразования широты и долготы в пиксели карты будут выполняться с помощью кода google, где используются формулы описывающие проекцию Меркатора:
Подробнее про проекции можно почитать здесь.
Начиная с версии Sql Server 2008, для работы с пространственными данными поддерживается тип данных GEOMETRY и GEOGRAPHY.
Существующие картографические сервисы от Yandex, Google или OpenStreetMap картографию как предоставляли в виде PNG картинок, фиксированного размера, обычно 256x256 точек. Хотя теперь существуют сервисы, где тайлы формируются с помощью таких технологий как SVG или CANVAS. Будем рассматривать растровые тайлы в формат PNG (picture network graphic). PNG формат поддерживает прозрачность (указывается в альфа канале), благодаря этому возможно наложение тайлов друг на друга без абсолютного перекрытия, при отображении нескольких слоёв.
Хранилище тайлов
Для каждого масштаба хранится определённый набор тайлов. Для масштаба 0-ого уровня – 1 тайл:
Для масштаба 1-ого уровня четыре тайла 2 * 2:
для масштаба n хранится 2n * 2n тайлов. Количество тайлов с увеличением номера масштаба увеличивается экспоненциально.
Тайлы хранятся в файловой системе Web-сервера и отправляются на машину клиента по http запросу, например:
someurl/layer{Z}/{X}/{Y}.png
где Z,X,Y соответственно масштаб, X позиция тайла, Y позиция тайла. Например, по следующим url доступен тайл с изображением Троицкого моста в Санкт-Петербурге:
b.tile.openstreetmap.org/15/19144/9524.png
В этом url:
15 – номер масштаба;
19144– X позиция тайла;
9524 – Y позиция тайла.
Естественно, что в разных системах формат URL запроса отличается. Вместо номеров тайлов и масштаба, для запроса тайлов может быть использован QUAD-ключ. Файлы тайлов могут быть отправлены клиенту сервером прямо из файловой системы или с помощью http обработчиков. Будем рассматривать вариант с X, Y, Z.
Этапы подготовки тайлов
- Формирование списка тайлов по геометрии каждого объекта;
- Генерация тайлов для каждого объекта;
- Объединение тайлов, для получения уникального набора;
- Сохранение в файловую систему.
Используемые функции
Для реализации задачи понадобится функция формирования геометрии тайла по X, Y позициям тайла и номеру масштаба. Геометрия тайла в нашем случае это прямоугольник покрывающий тайл с координатами углов, выраженными в широте и долготе. Формирования такой геометрии можно реализовать в SQL функции или в функции SQL CLR. Разница во времени выполнения SQL CLR функции и обычной SQL функции не заметна. Код SQL CLR функции реализован в классе Coords2PixelConversion в прилагаемых исходных кодах.
Следующая геометрия является контуром этого тайла, то есть проходит по его границам. Координаты вершин здесь – это долгота и широта.
'POLYGON ((30.322265625 59.955010262062061, 30.322265625 59.949509172252277, 30.333251953125 59.949509172252277, 30.333251953125 59.955010262062061, 30.322265625 59.955010262062061))'
tile.GetTileBounds(@level int, @x int, @y int)
CREATE FUNCTION [tile].[GetTileBounds] (@level int, @x int, @y int)
RETURNS geometry
AS
BEGIN
DECLARE @res GEOMETRY = NULL
IF @level IS NOT NULL AND @x IS NOT NULL AND @y IS NOT NULL
BEGIN
DECLARE @n1 FLOAT = PI() - 2.0 * PI() * @y / POWER(2.0, @level);
DECLARE @n2 FLOAT = PI() - 2.0 * PI() * (@y + 1) / POWER(2.0, @level);
DECLARE @top FLOAT = (180.0 / PI() * ATAN(0.5 * (EXP(@n1) - EXP(-@n1))));
DECLARE @bottom FLOAT = (180.0 / PI() * ATAN(0.5 * (EXP(@n2) - EXP(-@n2))));
DECLARE @tileWidth FLOAT = 360 / CONVERT(float, POWER(2, @level))
DECLARE @left FLOAT = @tileWidth * @x - 180,
@right FLOAT = @tileWidth * (@x + 1) - 180
SET @res = geometry::STPolyFromText('POLYGON (('
+ LTRIM(STR(@left, 25, 16)) + ' ' + LTRIM(STR(@top, 25, 16)) + ', '
+ LTRIM(STR(@left, 25, 16)) + ' ' + LTRIM(STR(@bottom, 25, 16)) + ', '
+ LTRIM(STR(@right, 25, 16)) + ' ' + LTRIM(STR(@bottom, 25, 16)) + ', '
+ LTRIM(STR(@right, 25, 16)) + ' ' + LTRIM(STR(@top, 25, 16)) + ', '
+ LTRIM(STR(@left, 25, 16)) + ' ' + LTRIM(STR(@top, 25, 16))
+ '))', 0)
END
RETURN @res
END
Как использовать эту функцию будет рассказано далее в этой статье.
Рассмотрим методы формирования списка тайлов. Для разных геометрий, можно подобрать различные подходы к формированию списка тайлов.
Способ 1:
Воспользуемся таким обстоятельством: Если на одном масштабе объект не пересекается с тайлом, то на масштабе с большим номером 4 тайла, которые перекрывают проверяемый тайл, так же не пересекаются с объектом. То есть при переходе к следующему масштабу выполняем проверку тайла, только если тайл предыдущего масштаба пересекается с геометрией объекта. Так исключаются лишние проверки, которые выполняются в способе 2.Способ 2:
На самом деле этот способ – сценарий худшего случая. Для каждого масштаба по каждому объекту определять поднабор тайлов функцией GEOMETRY::STEnvelope() и проверять пересечение тайла из этого поднабора с объектом. Этот способ менее эффективен, особенно для объектов с большой площадью или длиной, так как проверяется большее количество тайлов.Способ 3:
Для каждого объекта формировать геометрию тайловой сетки, по пересечению сетки с геометрией объекта получить набор точек. Для каждой, полученной точки определить два тайла и добавить в итоговый список тайлов.Например, для сложной географической линии проходящей через континент можно найти точки пересечения с сеткой проходящей по границам тайлов, и по этим точкам уже определить тайлы для рендеринга. Сетка создаётся в границах прямоугольной области, вмещающей линию, и представляет собой набор вертикальных и горизонтальных линий. Это гораздо эффективнее чем проверять каждый тайл в границах прямоугольной области объекта.
Опишем первый способ более подробно. Для геометрий объектов обладающих площадью набор тайлов для проверки пересечения с объектом может быть ограничен крайними тайлам прямоугольной области (bbox) перекрывающей объект.
По геометрии объекта (кроме типа геометрии POINT) формируется прямоугольник CLR функцией MSSQL GEOMETRY::STEnvelope(). Для объектов с геометрией POINT в качестве границы bbox используется прямоугольная область, перекрывающая иконку объекта на карте. Функция GetImageBound, возвращающая геометрию перекрывающую иконку реализована в классе GoogleProjection. Здесь же реализованы методы преобразования широты и долготы в номера позиций пикселей. Координаты угловых точек прямоугольной области выражены в широте и долготе. Далее мы получаем поднабор тайлов, покрывающих полученный прямоугольник. Для этого нам нужна функция преобразования географических координат в номер тайла на соответствующем масштабе. Для получения X и Y позиции тайла по долготе и широте можем использовать как SQL CLR функции, которые будут приведены далее, так и SQL функции приведённые ниже:
tile.GetXTilePos((@Longitude FLOAT, @Zoom INT)
tile.GetYTilePos((@Latitude FLOAT, @Zoom INT)
После определения позиций угловых тайлов, все тайлы, находящиеся в прямоугольной области между найденными угловыми тайлами проверяются на пересечение с геометрией объекта в функции tile.fn_FetchGeometryTilesZoomDepth().
CREATE FUNCTION tile.GetXTilePos(@Longitude FLOAT, @Zoom INT)
RETURNS INT
AS
BEGIN
DECLARE @D FLOAT,@E FLOAT,@F FLOAT,@G FLOAT, @tileY INT, @tileX INT
SET @D = 128 * POWER(2, @Zoom)
SET @E = ROUND(@D + @Longitude * 256 / 360 * POWER(2, @Zoom), 0)
SET @tileX = Floor(@E / 256);
RETURN @tileX
END
CREATE FUNCTION tile.GetYTilePos(@Latitude FLOAT, @Zoom INT)
RETURNS INT
AS
BEGIN
DECLARE @A FLOAT, @B FLOAT, @C FLOAT, @D FLOAT, @E FLOAT, @F FLOAT, @G FLOAT, @tileY INT
SET @D = 128 * POWER(2, @Zoom)
SET @A = Sin(PI() / 180 * @Latitude)
SET @B = -0.9999
SET @C = 0.9999
IF @A < @B SET @A = @B
IF @A > @C SET @A = @C
SET @F = @A
SET @G = Round(@D + 0.5 * Log((1.0 + @F) / (1.0 - @F)) * (-256) * POWER(2, @Zoom) / (2 * PI()),0)
SET @tileY = Floor(@G / 256)
RETURN @tileY
END
В функции tile.fn_FetchGeometryTilesZoomDepth() вычисляются левый верхний и правый нижний тайл минимальной прямоугольной области покрывающей геометрию. Затем для определенния пересечения фигуры с тайлом во вложенном цикле используем функцию tile.fn_GetTilesByTileNumZoomDepth() для каждого тайла в этой области проходя слева направо и сверху вниз от левого верхнего тайла до правого нижнего. Функция возвращает список тайлов, для которых было определено пересечение с геометрией объекта.
CREATE FUNCTION tile.fn_FetchGeometryTilesZoomDepth
( @GeoData GEOMETRY, @Zoom INT, @Depth INT)
RETURNS @retTiles TABLE ( Zoom INT, TileX INT, TileY INT)
AS
BEGIN
DECLARE @Left FLOAT, @Right FLOAT, @Top FLOAT, @Bottom FLOAT
DECLARE @CurrentXTile INT, @CurrentYTile INT, @Quanttiles INT
DECLARE @Envelope GEOMETRY, @RightTop GEOMETRY, @LeftBottom GEOMETRY
DECLARE @CurTileGeom GEOMETRY, @res GEOMETRY
DECLARE @tiletop FLOAT,@tilebottom FLOAT,@tileleft FLOAT, @tileright FLOAT
DECLARE @LeftTilePos INT,@RightTilePos INT,@TopTilePos INT
DECLARE @BottomTilePos INT
SET @envelope = @GeoData.STEnvelope()
SET @RightTop = @envelope.STPointN(3)
SET @LeftBottom = @envelope.STPointN(1)
SET @Right = @RightTop.STX
SET @Left = @LeftBottom.STX
SET @Top = @RightTop.STY
SET @Bottom = @LeftBottom.STY
SET @LeftTilePos = tile.GetXTilePos( @Left,@Zoom)
SET @RightTilePos = tile.GetXTilePos( @Right,@Zoom)
SET @TopTilePos = tile.GetYTilePos( @Top,@Zoom)
SET @BottomTilePos = tile.GetYTilePos( @Bottom,@Zoom)
SET @CurrentXTile = @LeftTilePos
WHILE @CurrentXTile <= @RightTilePos
BEGIN
SET @currentYTile = @TopTilePos
WHILE @CurrentYTile <= @BottomTilePos
BEGIN
INSERT INTO @retTiles (Zoom, TileX, TileY)
SELECT * FROM tile.fn_GetTilesByTileNumZoomDepth ( @GeoData, @Zoom, @CurrentXTile, @CurrentYTile, @Depth )
SET @CurrentYTile = @CurrentYTile + 1
END
SET @CurrentXTile =@CurrentXTile + 1
END
RETURN
END
Проверять пересечение геометрии тайла с геометрией объекта будем функцией GEOMETRY::STIntersects(). Если геометрия объекта и геометрия тайла пересекаются, то добавляем запись в предварительно созданную таблицу tile.TileOverlap и вызываем ту же функцию рекурсивно для четырёх тайлов следующего масштаба, покрывающие текущий, с параметром @Depth уменьшеным на еденицу. Проверка пересечения реализована в функции tile.fn_FetchGeometryTilesZoomDepth().
CREATE FUNCTION tile.fn_GetTilesByTileNumZoomDepth ( @GeoData GEOMETRY, @Zoom INT, @tileX INT, @tileY INT, @Depth INT)
RETURNS @retTiles TABLE ( Zoom INT, X INT, Y INT)
AS
BEGIN
DECLARE @currentTile TABLE ( Zoom INT, X INT, Y INT)
IF GEOGRAPHY::STGeomFromWKB([tile].[GetTileBounds](@Zoom, @tileX, @tileY).STAsBinary(),4326).STIntersects(GEOGRAPHY::STGeomFromWKB(@GeoData.MakeValid().STUnion(@GeoData.STStartPoint()).STAsBinary(),4326)) = 1
BEGIN
INSERT INTO @currentTile SELECT @Zoom , @tileX , @tileY
INSERT INTO @retTiles SELECT d.zoom, d.X, d.Y FROM @currentTile c
CROSS APPLY (SELECT * FROM [tile].[fn_GetTilesForObjectByTileNumZoomDepth]( @GeoData , c.Zoom + 1, c.X * 2, c.Y * 2, @Depth - 1) WHERE @Depth > 0) AS d
INSERT INTO @retTiles SELECT d.zoom, d.X, d.Y FROM @currentTile c
CROSS APPLY (SELECT * FROM [tile].[fn_GetTilesForObjectByTileNumZoomDepth]( @GeoData , c.Zoom + 1, c.X * 2 + 1, c.Y * 2, @Depth - 1) WHERE @Depth > 0) AS d
INSERT INTO @retTiles SELECT d.zoom, d.X, d.Y FROM @currentTile c
CROSS APPLY (SELECT * FROM [tile].[fn_GetTilesForObjectByTileNumZoomDepth]( @GeoData , c.Zoom + 1, c.X * 2, c.Y * 2 + 1, @Depth - 1) WHERE @Depth > 0) AS d
INSERT INTO @retTiles SELECT d.zoom, d.X, d.Y FROM @currentTile c
CROSS APPLY (SELECT * FROM [tile].[fn_GetTilesForObjectByTileNumZoomDepth]( @GeoData , c.Zoom + 1, c.X * 2 + 1, c.Y * 2 + 1, @Depth - 1) WHERE @Depth > 0) AS d
INSERT INTO @retTiles SELECT * FROM @currentTile
END
RETURN
END
Если необходимо сформировать тайлы для одного объекта, номера тайлов можно писать сразу в таблицу tile.Tile, так как набор тайлов будет уникальным. Для формирования тайла, с которым пересекаются геометрии нескольких объектов понадобиться объединение тайлов, созданных для разных объектов и накладывающихся друг на друга.
Функция tile.fn_GetTilesByTileNumZoomDepth() выполняет проверку пересечения геометрии объекта с тайлами масштабов, с учётом заданной глубины. Параметр @Depth указывает глубину проверки, если например @Zoom = 2 и @Depth = 1 то будет проверен тайл масштаба 2, и в случае наличия пересечения будут проверены 4 тайла масштаба 3. Необходимо проверить эти тайлы, так как они перекрывают тайл с предыдущего масштаба. Проверка пересечения должна выполнятся после преобразования типа данных GEOMETRY в GEOGRAPHY, это важно, так как для типа данных GEOGRAPHY проверка выполняется с учётом того, что все координаты точек геометрии в проекции 4326, то есть имеем дело с геометрическими объектами на сфере.
Пример с ломаной линией
Допустим, что хотим получить тайлы для ломанной соединяющей центр Санкт-Петербурга с центром Москвы. Протяжённость примерно 800 км. Ломанная будет проходить через населённые пункты: Новгород — Вышний волочок – Тверь.
'LINESTRING( 30.381113 59.971474, 31.26002 58.539215, 34.564158 57.591722, 35.915476 56.876838,37.622242 55.773125)'
Для этой геометрии с 3 по 17 масштаб получим всего тайлов 11076, распределение количества тайлов пересекающихся с геометрией по масштабам приведено в таблице 1
| Масштаб | Количество тайлов |
|---|---|
| 3 | 1 |
| 4 | 2 |
| 5 | 3 |
| 6 | 4 |
| 7 | 7 |
| 8 | 12 |
| 9 | 23 |
| 10 | 45 |
| 11 | 88 |
| 12 | 174 |
| 13 | 347 |
| 14 | 692 |
| 15 | 1383 |
| 16 | 2766 |
| 17 | 5529 |
Тайлы полученные для 3-ого и 4-ого масштаба показаны на рисунке 1:
Рисунок 1 – Тайлы: 3/4/2 и 4/9/4
Для каждого масштаба формируется поднабор тайлов и проверяется каждый без исключения тайл. На 4- 5 масштабе количество тайлов попадающих в прямоугольную область, полученную функцией GEOMETRY::STEnvelope() по геометрии объекта будет невелико. Всего тайлов на 4 масштабе 2^4*2^4 = 256. Но на 16 и 17 масштабе необходимо будет проверить на много больше тайлов. Исключение «лишних» проверок в первом способе ускорит работу. Для объектов с геометрией точка (POINT), оба способа будут иметь одинаковую эффективность.
Проверка пересечения
Функция GEOMETRY::STIntersects() может не определить пересечение геометрии объекта с геометрией тайла, так как функция STIntersects() для типа данных GEOMETRY работает с координатами на плоскости, а широта и долгота не являются декартовыми координатами. Поэтому, для достоверного определения пересечения конвертируем тип GEOMETRY в тип GEOGRAPHY. Следует обратить внимание, что в отличие от типа данных GEOMETRY тип данных GEOGRAPHY требует соблюдения направленности колец полигонов. Координаты внешних колец должны идти против часовой стрелки. Для внутренних колец (пустот) координаты должны перечисляться по часовой стрелке. Чтобы избежать ошибки при формировании географии используем фунции GEOMETRY::MakeValid() и GEOMETRY::STUnion() для получения корректной последовательности координат, при преобразования типа GEOMETRY в тип GEOGRAPHY. При создании географии указываем параметр SRID = 4326, это значит, что все операции над координатами выполняются в сферической системе.
На данном этапе, когда получен список тайлов, то есть таблицу с тремя колонками: Z, X, Y; дальнейшую работу можно выполнить с помощью программы для редеринга mapnik, которая позволяет выполнить формирование тайлов с изображением объектов. Организация доступа mapnik к базе данных Microsoft SQL Server требует определённых усилий. Подготовка к генерации тайлов в mapnik включает следующие шаги:
• Объявить стиль объектов для рендеринга;
• Описать источник данных геометрии объектов;
• Указать таблицу со списком тайлов для рендеринга, чтобы не генерировать все тайлы подряд.
Мы будем выполнять генерацию тайлов внутри базы данных MS SQL Server 2008. Для этого необходимо реализовать несколько CLR функций для работы с геометрией, хранящейся в базе данных. Основные функции, которые нам понадобятся приведены ниже:
- tile.IconTile(),
- tile.ShapeTile(),
- tile.TileAgg(),
- tile.SaveToFolderByZoomXY()
.
Таблицы для хранения образов тайлов
На рисунке 2 показана структура таблиц, где сохраняется список тайлов для рендеринга. В поле Data в этих таблицах будет сохраняться PNG картинка с изображением объектов попадающих на тайл. Хранение и обработка большого количества картинок в таблице может сказаться на производительности. Для данной задачи более подходящим вариантом является формирование образов тайлов вне базы данных по сформированному в таблице списку тайлов по каждому объекту и последующее сохранение в файловую систему.
Рисунок 2 – Таблицы для хранения списка тайлов
Размещение иконки на тайле
Разберём алгоритмы позиционирования иконки на тайле (тип геометрии POINT).
Есть широта и долгота некоторого объекта, есть список тайлов текущего масштаба, на которые накладывается иконка. Формирование перечня тайлов было описано ранее. Вычисление позиции иконки на тайле состоит из следующих действий:
1. Сначала преобразуем широту и долготу в абсолютные пиксельные координаты;
2. Затем, для каждого тайла, из имеющегося списка, на текущем масштабе вычисляем абсолютные пиксельные координаты левого верхнего угла. координаты левого верхнего пикселя тайла (pixXtile,pixYtile) вычисляем умножением номера x и y позиции тайла на его размер, в нашем случае это 256 пикселей;
3. Разность между абсолютными пиксельными координатами объекта и абсолютными пиксельными координатами левого верхнего угла тайла определяются в функция GetPixelXOnTile() и GetPixelXOnTile(). Эта разность есть относительные пиксельные координаты центра иконки на тайле;
4. Чтобы отрисовать иконку на тайле нужно получить границы области отрисовки на тайле в пикселях, в которую будет происходить вставка. Относительные пиксельные координаты объекта на тайле получены на предыдущем шаге. Теперь по размеру иконки определяются границы прямоугольной области для вставки.
5. Выполняем отрисовку иконки на тайле.
[Microsoft.SqlServer.Server.SqlFunction]
public static SqlBinary IconTile(SqlBinary image, SqlInt32 zoom, SqlDouble Lon, SqlDouble Lat, SqlInt32 xTile, SqlInt32 yTile, SqlDouble scale)
{
SqlBinary result = null;
using (Icon2TileRendering paster = new Icon2TileRendering())
{
using (MemoryStream ms = new MemoryStream())
{
ms.Write(image.Value, 0, image.Length);
SetBeginPosition(ms);
paster.PasteFromStreamScaledImageToTile((int)zoom, (double)Lon, (double)Lat, (int)xTile, (int)yTile, (double)scale, ms);
result = paster.GetBytes();
}
}
return result;
}
#region [Pixel Position Calculation]
Rectangle GetTargetBound(int zoom, double Lon, double Lat, int xTile, int yTile, int width, int height)
{
int xPix = _conv.FromLongitudeToXPixel(Lon, zoom);
int yPix = _conv.FromLatitudeToYPixel(Lat, zoom);
int xPos = xPix - xTile * TILE_SIZE;
int yPos = yPix - yTile * TILE_SIZE;
int halfWidth = width / 2;
int halfHeight = height / 2;
return new Rectangle(xPos - halfWidth, yPos - halfHeight, width, height
}
int GetPixelXOnTile(int zoom, double Lon, int xTile)
{
return _conv.FromLongitudeToXPixel(Lon, zoom) - xTile * TILE_SIZE;
}
int GetPixelYOnTile(int zoom, double Lat, int yTile)
{
return _conv.FromLatitudeToYPixel(Lat, zoom) - yTile * TILE_SIZE;
}
#endregion [Pixel Position Calculation]
/// <summary>
/// Размещает иконку на тайле
/// </summary>
/// <param name="zoom"></param>
/// <param name="Lon"></param>
/// <param name="Lat"></param>
/// <param name="xTile"></param>
/// <param name="yTile"></param>
/// <param name="iconImage"></param>
public void PasteImageToTileByLatLon(int zoom, double Lon, double Lat, int xTile, int yTile, Bitmap iconImage)
{
int width = iconImage.Width;
int height = iconImage.Height;
CopyRegionIntoImage(iconImage, new Rectangle(0, 0, width, height), GetTargetBound(zoom, Lon, Lat, xTile, yTile, width, height));
}
Объединение тайлов
На один и тот же тайл могут накладываться иконки нескольких объектов. Чтобы получить тайлы со всеми объектами, можно сначала создать тайлы для каждого объекта, затем объединить их в один. Такое решение можно реализовать с помощью группировки строк таблицы базы данных, для этого создана CLR функция агрегат tile.TileAgg(), объединяющая тайлы в один. Это решение имеет один минус, так как на каждый объект пресекающийся с тайлом мы будем хранить отдельную запись с BINARY полем, хранящим картинку тайла, с избражением только данного объекта, что требует большого количества памяти. Более правильное решение – использовать один экземпляр тайла, и последовательно отображать на нем все попадающие на него иконки объектов. Таким образом, расходуем меньше памяти. В этом случае группироать просто нечего. Мы хотим использовать группировку.
CREATE PROC [tile].[FillObjectTilesIntersection]( @StartZoom INT, @EndZoom INT)
AS
BEGIN
DECLARE @CurrentZoom INT
SET @CurrentZoom = @StartZoom
WHILE @CurrentZoom <= @EndZoom
BEGIN
INSERT INTO tile.Tile (X,Y,Data,Zoom)
SELECT t.TileX,t.TileY, [tile].[TileAgg]
(tile.IconTile(i.Data, @CurrentZoom,o.Longitude,o.Latitude,t.tileX,t.tileY, 1.0)
),@CurrentZoom AS Zoom
FROM tile.Shape o
INNER JOIN tile.[Image] i ON i.ImageID = o.ImageID
CROSS APPLY tile.fn_FetchObjectTiles(tile.GetIconBoundForZoom(o.Longitude, o.Latitude, 64, 64, @CurrentZoom, 0),@CurrentZoom) t
WHERE o.TypeID = @TypeID
GROUP BY t.TileX,t.TileY
SET @CurrentZoom = @CurrentZoom + 1
END
END
В качестве источника объектов будем использовать таблицу tile.Object с координатами объектов и идентификатором образа иконки, хранящейся в таблице tile.Image в поле типа Binary.
DECLARE @Image VARBINARY(MAX)
SELECT TOP (1) @image = ( SELECT * FROM OPENROWSET(BULK N'd:\media\icons\pin_orange.png', SINGLE_BLOB) As tile)
SELECT [tile].[SaveToFolderByZoomXY]([tile].[IconTile](@image, 3,30.381113, 59.971474, 4, 2, 1.0), N'D:\Tiles\',3,4,2)
SELECT [tile].[SaveToFolderByZoomXY]([tile].[IconTile](@image, 4,30.381113, 59.971474, 9, 4, 1.0), N'D:\Tiles\',4,9,4)
Отрисовка геометрии на тайле
Для ломаной линии (POLYLINE, MULTIPOLYLINE) объединяем геометрию тайла с геометрией полилинии, таким образом часть ломаной находящейся вне области тайла исключается. Алгоритм определения контура и закрашиваемой области, может быть применён для геометрий обладающих площадью, то есть POLYGON, MULTIPOLYGON, GEOMETRYCOLLECTION содержащая POLYGON или MULTYPOLYGON. Алгоритм реализован в классе ShapeToTileRendering и включает следующие этапы:
1. Координаты (широта, долгота) геометрии преобразуются в пиксельные координаты с учётом масштаба по формулам преобразования широты, долготы в пиксели PSG3857 (Google проекция) и вычитаем из каждой координаты полученной геометрии координаты левого верхнего пикселя целевого тайла. Получаем, так называемую геометрию (A). Эти действия реализованы в функции ConvertToZoomedPixelZeroedByTileGeometry(poly, Z, X, Y)
2. Формируется геометрия (B) тайла в пиксельных координатах с учётом масштаба
3. Формируется геометрия (C) полученная пересечением (STIntersection) пиксельной геометрии тайла (B) с геометрией объекта (A)
4. Формируется геометрия (D) полученная в результате пересечения контура геометрии (C) и контура геометрии тайла (B), получаем линии проходящие по границе тайла и граничащие с предполагаемой закрашиваемой областью полигона внутри тайла. Формируется геометрия (E) полученная в результате вычитания, с помощью функции .STDifference ( other_geometry )
5. Геометрия (E) и есть контур для отрисовки, который получается вычитанием из контура(LINESTRING или MULTILINSTRING) геометрии (С) геометрии (D) с помощью функции
6. Заливается полигон геометрии (С) – получена заливаемая область
7. Рисуется геометрия (E) -контур полигона после исключения пересечения с границами тайла
8. Повторяем шаги с 1 по 7 для всех тайлов текущего объекта и сохраняем их в таблицу tile.TileOverlap
Рассмотрим первые 3 этапа более подробно на примере тайла 15-ого масштаба в X позиции с номером 19144 и Y позиции с номером 9524. Выполним вышеперечисленные операции с помощью T-SQL скриптов. Для начала получаем геометрию границ тайла, использую следующий скрипт:
SELECT [tile].[GetTileBounds](15,19144,9524).ToString()
Получили следующий результат:
'POLYGON ((30.322265625 59.955010262062061, 30.322265625 59.949509172252277, 30.333251953125 59.949509172252277, 30.333251953125 59.955010262062061, 30.322265625 59.955010262062061))'
Далее построим геометрию в виде ромба например, которая будет пересекаться с выбранным тайлом. Пусть центры геометрии и тайла будут совпадать. Для этого воспользуемся функцией построения географического сектора на условной поверхности земли. Функция реализована, исходя из допущения, что земля имеет форму сферы с радиусом 6367 км. Параметрами функции являются долгота и широта центра дуги, азимут(направление сектора), угол разворота дуги, радиус в метрах и шаг приращения угла в градусах. Чем меньше приращение, тем более точная дуга и больше точек в геометрии. Если указать нулевой азимут, угол разворота равным 360 градусов и приращение 90 градусов, то получим наш ромб. В функции используется цикл, где c помощью формул сферической тригонометрии выполняются приращения к углу и определяются координаты точек на дуге окружности, лежащей на поверхности земли. Функция возвращает геометрию полигона на основе этих точек:
SELECT [tile].[fn_GetCircleSegment](30.3277587890625, 59.952259717159905,0,360,440,90)
CREATE FUNCTION [tile].[fn_ GetCircleSegment]
(@X float, @Y float, @azimuth float, @angle float, @distance float, @step FLOAT)
RETURNS geometry
WITH EXEC AS CALLER
AS
BEGIN
IF @X IS NULL OR @Y IS NULL OR @azimuth IS NULL OR ISNULL(@angle, 0) = 0 OR ISNULL(@distance, 0) <= 0
RETURN NULL
DECLARE @sectorStepAngle FLOAT
SET @sectorStepAngle = @step
IF ABS(@angle) > 360
SET @angle = 360
DECLARE @pointsStr VARCHAR(MAX)
DECLARE @firstPointsStr VARCHAR(MAX)
DECLARE @earthRadius FLOAT
DECLARE @lat FLOAT
DECLARE @lon FLOAT
DECLARE @d FLOAT
IF ABS(@angle) < 360
SET @pointsStr = LTRIM(STR(@X, 25, 16)) + ' ' + LTRIM(STR(@Y, 25, 16))
ELSE SET @pointsStr = ''
SET @earthRadius = 6367
SET @lat = RADIANS(@Y)
SET @lon = RADIANS(@X)
SET @d = (@distance / 1000) / @earthRadius
DECLARE @angleStart FLOAT
DECLARE @angleEnd FLOAT
SET @angleStart = @azimuth - @angle / 2;
SET @angleEnd = @azimuth + @angle / 2;
DECLARE @pointsCount INT
SET @pointsCount = FLOOR(@angle / @sectorStepAngle)
DECLARE @brng FLOAT
DECLARE @latRadians FLOAT
DECLARE @lngRadians FLOAT
DECLARE @ptX FLOAT
DECLARE @ptY FLOAT
DECLARE @i INT
SET @i = 0
DECLARE @addPrefix TINYINT
IF ABS(@angle) < 360
SET @addPrefix = 1
ELSE SET @addPrefix = 0
WHILE @i <= @pointsCount
BEGIN
SET @brng = RADIANS(@angleStart + @i * @sectorStepAngle);
SET @latRadians = ASIN(SIN(@lat) * COS(@d) + COS(@lat) * SIN(@d) * COS(@brng));
SET @lngRadians = @lon + ATN2(SIN(@brng) * SIN(@d) * COS(@lat), COS(@d) - SIN(@lat) * SIN(@latRadians));
SET @ptX = 180.0 * @lngRadians / PI();
SET @ptY = 180.0 * @latRadians / PI();
IF @addPrefix = 1
BEGIN
SET @pointsStr += ', ' + LTRIM(STR(@ptX, 25, 16)) + ' ' + LTRIM(STR(@ptY, 25, 16))
END
ELSE
BEGIN
SET @pointsStr += LTRIM(STR(@ptX, 25, 16)) + ' ' + LTRIM(STR(@ptY, 25, 16))
SET @addPrefix = 1
END
IF @i = 0
SET @firstPointsStr = LTRIM(STR(@ptX, 25, 16)) + ' ' + LTRIM(STR(@ptY, 25, 16))
IF @i = @pointsCount AND (@angleStart + @pointsCount * @sectorStepAngle) < @angleEnd
BEGIN
SET @brng = RADIANS(@angleEnd);
SET @latRadians = ASIN(SIN(@lat) * COS(@d) + COS(@lat) * SIN(@d) * COS(@brng));
SET @lngRadians = @lon + ATN2(SIN(@brng) * SIN(@d) * COS(@lat), COS(@d) - SIN(@lat) * SIN(@latRadians));
SET @ptX = 180.0 * @lngRadians / PI();
SET @ptY = 180.0 * @latRadians / PI();
SET @pointsStr = @pointsStr + ', ' + LTRIM(STR(@ptX, 25, 16)) + ' ' + LTRIM(STR(@ptY, 25, 16))
END
SET @i = @i + 1
END
IF ABS(@angle) < 360
SET @pointsStr += ', ' + LTRIM(STR(@X, 25, 16)) + ' ' + LTRIM(STR(@Y, 25, 16))
ELSE
SET @pointsStr += ', ' + @firstPointsStr
RETURN geometry::STPolyFromText('POLYGON ((' + @pointsStr + '))', 0)
END
GO
В качестве альтернативы, написана CLR функция, выполняющая те же вычисления. Можно заметить что значительного отличия во времени выполнения этих функций не наблюдается.
/// <summary>
/// Формирует геометрию сектора
/// </summary>
/// <param name="longitude"></param>
/// <param name="latitude"></param>
/// <param name="azimuth"></param>
/// <param name="angle"></param>
/// <param name="radius"></param>
/// <returns></returns>
[Microsoft.SqlServer.Server.SqlFunction]
public static SqlGeometry DrawGeoSpatialSectorVarAngle(SqlDouble longitude, SqlDouble latitude, SqlDouble azimuth,
SqlDouble angle, SqlDouble radius, SqlDouble stepAngle)
{
if (longitude == SqlDouble.Null || latitude == SqlDouble.Null || azimuth == SqlDouble.Null ||
angle == SqlDouble.Null || radius == SqlDouble.Null || radius == 0 || angle == 0)
return SqlGeometry.Parse("GEOMETRYCOLLECTION EMPTY");
SqlGeometryBuilder builder = new SqlGeometryBuilder();
builder.SetSrid(0);
builder.BeginGeometry(OpenGisGeometryType.Polygon);
double firstPointLon;
double firstPointLat;
double sectorStepAngle = (double) stepAngle;
const double earthRadius = 6367.0;
double lat = (double) latitude;
double lon = (double) longitude;
double azim = (double) azimuth;
double ang = (double) angle;
double piRad = (Math.PI/180.0);
double tLat = piRad*lat;
double tLon = piRad*lon;
double distkm = ((double) radius/1000)/earthRadius;
double angleStart = azim - ang/2;
double angleEnd = azim + ang/2;
var _angle = Math.Abs(ang);
if (_angle > 360.0)
{
angle = 360.0;
}
int pointCount = (int) Math.Floor(ang/sectorStepAngle);
double brng;
double latRadians;
double lngRadians;
double ptX;
double ptY;
int i = 0;
if (angle < 360.0)
{
builder.BeginFigure(lon, lat);
firstPointLon = lon;
firstPointLat = lat;
}
else
{
brng = piRad*(angleStart);
latRadians = Math.Asin(Math.Sin(tLat)*Math.Cos(distkm) + Math.Cos(tLat)*Math.Sin(distkm)*Math.Cos(brng));
lngRadians = tLon + Math.Atan2(Math.Sin(brng)*Math.Sin(distkm)*Math.Cos(tLat),
Math.Cos(distkm) - Math.Sin(tLat)*Math.Sin(latRadians));
ptX = 180.0*lngRadians/Math.PI;
ptY = 180.0*latRadians/Math.PI;
builder.BeginFigure(ptX, ptY);
firstPointLon = ptX;
firstPointLat = ptY;
}
while (i <= pointCount)
{
brng = piRad*(angleStart + i*sectorStepAngle);
latRadians = Math.Asin(Math.Sin(tLat)*Math.Cos(distkm) + Math.Cos(tLat)*Math.Sin(distkm)*Math.Cos(brng));
lngRadians = tLon + Math.Atan2(Math.Sin(brng)*Math.Sin(distkm)*Math.Cos(tLat),
Math.Cos(distkm) - Math.Sin(tLat)*Math.Sin(latRadians));
ptX = 180.0*lngRadians/Math.PI;
ptY = 180.0*latRadians/Math.PI;
builder.AddLine(ptX, ptY);
i = i + 1;
}
if (((angleStart + pointCount * sectorStepAngle) < angleEnd))
{
brng = piRad * (angleEnd);
latRadians = Math.Asin(Math.Sin(tLat) * Math.Cos(distkm) + Math.Cos(tLat) * Math.Sin(distkm) * Math.Cos(brng));
lngRadians = tLon + Math.Atan2(Math.Sin(brng) * Math.Sin(distkm) * Math.Cos(tLat),
Math.Cos(distkm) - Math.Sin(tLat) * Math.Sin(latRadians));
ptX = 180.0 * lngRadians / Math.PI;
ptY = 180.0 * latRadians / Math.PI;
builder.AddLine(ptX, ptY);
}
builder.AddLine(firstPointLon, firstPointLat);
builder.EndFigure();
builder.EndGeometry();
return builder.ConstructedGeometry;
}
DECLARE @bbox GEOMETRY
DECLARE @octagon GEOMETRY
SELECT @bbox = [tile].[GetTileBounds](15,19144,9524),
@octagon = [tile].[fn_GetCircleSegment](30.3277587890625, 59.952259717159905,0,360,440,90)
Здесь 30.3277587890625, 59.952259717159905 – координаты центра тайла;
Получим геометрию пересечения, пока ещё с широтой и долготой в координатах:
SELECT @bbox.STIntersection(@octagon)
Получаем следующую геометрию:
'POLYGON ((30.3253442162734 59.949509172234684,
30.3301733618516 59.949509172234684,
30.333251953125 59.9510505967796,
30.333251953125 59.953468509045528,
30.330173073498937 59.955010262085125,
30.325344504626063 59.955010262085125,
30.322265625 59.953468509045528,
30.322265625 59.9510505967796,
30.3253442162734 59.949509172234684))'
Будем использовать функции преобразования широты и долготы в пиксели, чтобы получить порядковые номера пикселей на виртуальном плоском полотне карты:
SELECT [tile].[GetPixelXPosFromLongitude](30.3253442162734,15), [tile].[GetPixelYPosFromLatitude](59.949509172234684,15)
, [tile].[GetPixelXPosFromLongitude](30.3301733618516,15), [tile].[GetPixelYPosFromLatitude]( 59.949509172234684,15)
, [tile].[GetPixelXPosFromLongitude](30.333251953125,15), [tile].[GetPixelYPosFromLatitude]( 59.9510505967796,15)
, [tile].[GetPixelXPosFromLongitude](30.333251953125,15), [tile].[GetPixelYPosFromLatitude]( 59.953468509045528,15)
, [tile].[GetPixelXPosFromLongitude](30.330173073498937,15), [tile].[GetPixelYPosFromLatitude]( 59.955010262085125,15)
, [tile].[GetPixelXPosFromLongitude](30.325344504626063,15), [tile].[GetPixelYPosFromLatitude]( 59.955010262085125,15)
,[tile].[GetPixelXPosFromLongitude](30.322265625,15), [tile].[GetPixelYPosFromLatitude]( 59.953468509045528,15)
, [tile].[GetPixelXPosFromLongitude](30.322265625,15), [tile].[GetPixelYPosFromLatitude]( 59.9510505967796,15)
, [tile].[GetPixelXPosFromLongitude](30.3253442162734,15), [tile].[GetPixelYPosFromLatitude]( 59.949509172234684,15)
| Долгота | Широта | X пиксель 15-ого масштаба | Y пиксель 15-ого масштаба |
|---|---|---|---|
| 30.3253442162734 | 59.949509172234684 | 4900936 | 2438400 |
| 30.3301733618516 | 59.949509172234684 | 4901048 | 2438400 |
| 30.333251953125 | 59.9510505967796 | 4901120 | 2438328 |
| 30.333251953125 | 59.953468509045528 | 4901120 | 2438216 |
| 30.330173073498937 | 59.955010262085125 | 4901048 | 2438144 |
| 30.325344504626063 | 59.955010262085125 | 4900936 | 2438144 |
| 30.322265625 | 59.953468509045528 | 4900864 | 2438216 |
| 30.322265625 | 59.9510505967796 | 4900864 | 2438328 |
| 30.3253442162734 | 59.949509172234684 | 4900936 | 2438400 |
По полученным пиксельным координатам сформируем геометрию контура пересечения тайла и геометрии объекта:
SELECT GEOMETRY::STGeomFromText('LINESTRING(4900936 2438400, 4901048 2438400, 4901120 2438328, 4901120 2438216, 4901048 2438144, 4900936 2438144, 4900864 2438216, 4900864 2438328, 4900936 2438400
)',0)
Геометрия (С) формируемая на шаге 3, область подлежащая заливке, выделена зелёным на рисунке 4 ниже.
Рисунок 4 – общая область геометрии тайла и геометрии объекта
Геометрия (D) нам не особенно нужна, больше нас интересует набор линий контура для отрисовки на тайле, его получаем, как описано в шаге 5, вычитанием из контура геометрии (С) геометрии (D), то что получилось выделено синим на рисунке 5.
Таким образом мы получили область для заливки и набор линий контура для данного тайла.
Следующая геометрия есть геометрия (E):
SELECT GEOMETRY::STGeomFromText('MULTILINESTRING((4901048 2438400, 4901120 2438328),( 4901120 2438216, 4901048 2438144),( 4900936 2438144, 4900864 2438216), (4900864 2438328, 4900936 2438400)
)',0)
Рисунок 5 – Контур геометрии для рендеринга (выделен синим)
Следующий T-SQL скрипт создаёт изображения тайлов в формате PNG и сохраняет в файловую систему по папкам соответствующим дереву Z/X/Y. При отсутствии папки создаются.
DECLARE @bbox GEOMETRY
DECLARE @rhomb GEOMETRY
DECLARE @image VARBINARY(MAX)
SELECT @bbox = [tile].[GetTileBounds](15,19144,9524), @rhomb = [tile].[fn_GetSector](30.3277587890625, 59.952259717159905,0,360,440,90)
SET @image = [tile].[ShapeTile]( @octagon,15,19144,9524,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19144,9524)
SET @image = [tile].[ShapeTile]( @octagon,15,19143,9524,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19143,9524)
SET @image = [tile].[ShapeTile]( @octagon,15,19145,9524,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19145,9524)
SET @image = [tile].[ShapeTile]( @octagon,15,19144,9523,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19144,9523)
SET @image = [tile].[ShapeTile]( @octagon,15,19144,9525,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19144,9525)
Полученные PNG файлы на рисунке ниже:
Для создание одного тайла вызываем метод DrawPartObjectShapeOnTile() класса ShapeToTileRendering:
/// <summary>
/// Отрисовка части геометриии на тайле
/// </summary>
/// <param name="shape">Полная геометрия объекта с градусами в качестве координат</param>
/// <param name="X">X позиция тайла</param>
/// <param name="Y">Y позиция тайла</param>
/// <param name="Zoom">номер уровня масштаба</param>
/// <param name="argbFill">цвет заполнения в формате ARGB</param>
/// <param name="argbStroke">цвет контура</param>
/// <param name="strokeWidth">ширина контура</param>
public void DrawPartObjectShapeOnTile(SqlGeometry shape, int X, int Y, int Zoom, string argbFill,
string argbStroke, int strokeWidth)
{
PasteShapeOnTile(CreateColor(argbFill), CreateColor(argbStroke), strokeWidth,
CutPolygonByZoomedPixelZeroTile(shape, X, Y, Zoom));
}
В метотде PasteSha
