Рисуем на тайлах электронной карты в MSSQL

Хочу рассказать читателям хабра-сообщества как используя CLR библиотеку Microsoft.SqlServer.Types можно формировать тайлы для электронной карты. В статье пойдёт речь о генерации списка картографических тайлов для их дальнейшего рендеринга. Будет описан алгоритм генерации тайлов по геометрии объектов, хранящейся в базе данных MS SQL 2008. Весь процесс рендеринга шаг за шагом будет рассмотрен на примере в конце статьи.

6a3a686e21ca4a818baa1248891b6e6c.png

Содержание



Проблема
Исходные данные
Решение
Хранилище тайлов
Этапы подготовки тайлов
Используемые функции
Пример с ломаной линией
Проверка пересечения
Таблицы для хранения образов тайлов
Размещение иконки на тайле
Объединение тайлов
Отрисовка геометрии на тайле
Заключение

Проблема



Когда в браузере отображается большое количество гео данных в векторной графике (с помощью 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 тайл:
image

Для масштаба 1-ого уровня четыре тайла 2 * 2:
a192980365ba4a999d4260035d96d4ed.jpg
для масштаба 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))'



Код скалярной SQL функции tile.GetTileBounds()
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().

SQL Функция получения X позиции тайла для долготы и номера масштаба
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


Функция получения Y позиции тайла для широты и номера масштаба
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

Таблица 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:
c84948b5e595451caf72736ee6995ae3.pngc6afd548047b418a970c57c2ff2ae7d0.png

Рисунок 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 картинка с изображением объектов попадающих на тайл. Хранение и обработка большого количества картинок в таблице может сказаться на производительности. Для данной задачи более подходящим вариантом является формирование образов тайлов вне базы данных по сформированному в таблице списку тайлов по каждому объекту и последующее сохранение в файловую систему.
010700c06a0741af99006ddb3310e58b.png
Рисунок 2 – Таблицы для хранения списка тайлов

Размещение иконки на тайле



Разберём алгоритмы позиционирования иконки на тайле (тип геометрии POINT).
Есть широта и долгота некоторого объекта, есть список тайлов текущего масштаба, на которые накладывается иконка. Формирование перечня тайлов было описано ранее. Вычисление позиции иконки на тайле состоит из следующих действий:
1. Сначала преобразуем широту и долготу в абсолютные пиксельные координаты;
2. Затем, для каждого тайла, из имеющегося списка, на текущем масштабе вычисляем абсолютные пиксельные координаты левого верхнего угла. координаты левого верхнего пикселя тайла (pixXtile,pixYtile) вычисляем умножением номера x и y позиции тайла на его размер, в нашем случае это 256 пикселей;
3. Разность между абсолютными пиксельными координатами объекта и абсолютными пиксельными координатами левого верхнего угла тайла определяются в функция GetPixelXOnTile() и GetPixelXOnTile(). Эта разность есть относительные пиксельные координаты центра иконки на тайле;
4. Чтобы отрисовать иконку на тайле нужно получить границы области отрисовки на тайле в пикселях, в которую будет происходить вставка. Относительные пиксельные координаты объекта на тайле получены на предыдущем шаге. Теперь по размеру иконки определяются границы прямоугольной области для вставки.
5. Выполняем отрисовку иконки на тайле.

CLR SQL функция размещающая иконку на тайле
[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.

Cкрипт формирования тайлов 3/4/2 и 4/9/4 с позиционированием иконки на долготе 30.381113 и широте 59.971474
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)



Рисунок 3 – Полученные тайлы с иконкой
52ce9758bcc84a00b61a8f80d8325966.pngddf973ad106e4baeaf45a45f166b6800.png

Отрисовка геометрии на тайле



Для ломаной линии (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 функция, выполняющая те же вычисления. Можно заметить что значительного отличия во времени выполнения этих функций не наблюдается.

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))'


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

Преобразование широты и долготы в X и Y позиции пикселей
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 ниже.
bfe63b6f342d45e7bdd3f1e29eb557c4.jpg
Рисунок 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)


3bdaa6d991df4709bc60306b313e7769.jpg
Рисунок 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 файлы на рисунке ниже:
c1eda1b7c5cd4e4590f9d4219e2b05e7.png
1612e81fd13e44eebd647b85569dbc82.png81badff4353b4cecac28cd7c0eb10cc2.pngc9530eb433244bbb8fc95d84c2acde94.pngbc07bca024c44606b56f87af40de525d.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

© Habrahabr.ru