目录

引言

初始数据和条件

解决方案

切片存储

准备切片的步骤

要使用的函数

折线对象形状示例

交集检测

用于存储切片图像的表

在特定位置的切片上绘制图标

切片组合

在切片上绘制几何形状

结论

引言

本文介绍如何在地理切片上绘制几何图形。当大量图标或形状以矢量格式（例如 SVG 或使用 CANVAS）在浏览器中的电子地图上显示时，会消耗客户端机器的大量资源。因此，我们将它们显示为栅格图像。我们将了解如何填充要渲染的切片列表。本文探讨了在切片上定位图标以及在切片上绘制图形基元的方法。它假设图形基元存储在 MS SQL 数据库表的 GEOMETRY 字段中。本文最后将通过一个实际示例逐步介绍对象形状的切片渲染整个过程。

初始数据和条件

几何形状的子集存储在 Microsoft SQL 2008 数据库中。顶点坐标为经纬度 [EPSG:4326]。空间数据字段的数据类型为 GEOMETRY。只有一个点作为形状的对象（具有 Point 形状类型的几何图形）在地图上显示为 PNG 格式的图标。如果对象具有一系列截距作为其形状，则将其绘制为具有指定宽度的折线。在数据库中，它具有 Linestring 或 Multilinestring 形状类型。多边形 GEOMETRY 在切片上绘制为带描边轮廓的填充形状。多面 GEOMETRY 是一组带指定宽度描边轮廓的填充形状。我们的切片系统应适应 Web 墨卡托投影 (http://spatialreference.org/ref/sr-org/7483/)

解决方案

让我们为数据库中的所有地理对象准备切片。这意味着对象将通过切片层显示，而不是在客户端上将其绘制为 PNG 图标或不同于 Point 的形状类型的矢量集。为了实现这一点，我们应该为所有对象的每个缩放级别创建切片。我们将使用 Web Mercator 投影和来自 GoogleMapsAPIProjection.cs 的 google-code。使用 Mercator 投影的公式，我们将经纬度转换为特定缩放级别地图上的切片像素编号。我们还使用这些公式计算切片位置，以选择与对象形状相交的所有切片。我们可以使用几种投影，要了解更多信息，请访问此链接：Mercator_projection。在 Microsoft Sql Server 2008 中引入了空间数据类型，即 GEOMETRY 和 GEOGRAPHY。有 OpenStreetMap、谷歌地图等地理服务。这些服务提供制图学作为固定大小的 PNG 图片，以便在 Web 浏览器中的地图上平铺。通常，该图片的正常大小为 256 像素。尽管现在一些地理服务提供 SVG 或 CANVAS 技术的切片。让我们使用“老式”PNG 格式的切片。PNG 支持透明度，每个像素都在 Alpha 通道中指定。它允许我们将切片堆叠在一起以在地图上显示。

切片存储

每个缩放级别都有固定数量的切片。对于零缩放级别，只有一个切片，对于缩放级别一，已有四个切片。对于缩放级别 n，我们有 2n * 2n 个切片。随着缩放级别编号的增加，切片的数量呈指数增长。下图显示了缩放级别一的切片

图 1 – 墨卡托投影缩放级别一的切片

切片根据其水平、垂直位置和缩放级别进行索引或编号，然后切片存储在 Web 服务器的文件系统中，并通过 http 发送到客户端机器。请求 URL 可能如下例所示

http://someurl/layer/{Z}/{X}/{Y}.png

其中 Z、X 和 Y 分别是缩放级别、切片的水平编号和垂直编号。例如，通过以下 URL，您可以下载俄罗斯圣彼得堡特罗伊茨基桥附近的切片

http://b.tile.openstreetmap.org/15/19144/9524.png

上述 URL 包含

15 - 缩放级别编号；

19144 切片的水平位置；

9524 – 切片的垂直位置。

切片请求 URL 的格式因系统而异。除了使用缩放级别编号和切片位置外，还可以使用 QUAD 键 (http://wiki.openstreetmap.org/wiki/QuadTiles)。切片文件通过服务器文件系统或 http 处理程序发送到客户端。让我们观察 Z、X、Y 的情况，其中在每个 X 列文件夹中，所有 Y 切片都以 PNG 格式存储。还假设在每个缩放级别文件夹中都有 X 列文件夹。因此，在顶层我们有缩放级别文件夹。在第二层有 X 文件夹，其名称与 X 切片位置匹配，第三层（最后一层）包含该 X 切片列的所有切片图片。切片文件的名称是切片的垂直位置。

准备切片的步骤

• 根据对象几何图形相交的切片编号形成列表；
• 根据上一步形成的列表为每个对象生成切片；
• 组合切片以获得唯一列表；
• 根据唯一的切片编号列表将切片保存到文件系统。

切片列表包含 Z、X、Y 编号和对象 ID。

要使用的函数

为了完成渲染任务，我们将实现一些函数。我们将编写的第一个函数是切片边界几何图形创建函数。让我们定义切片边界几何图形。它是一个矩形，其坐标可以是像素或经纬度。这里我们需要它的经纬度。我们的函数使用 X、Y 切片位置和缩放级别编号作为参数。函数返回一个 GEOMETRY 作为结果。结果是一个矩形，覆盖切片，并将其左上顶点和右下顶点作为其在经纬度中的坐标。换句话说，矩形的线沿着切片边界。函数可以写成普通的标量 SQL 函数或 SQL CLR 函数。它假设两者在执行时性能或处理器利用率没有区别。SQL CLR 函数在附件源代码的 Coords2PixelConversion 类中实现。

以下几何表示是缩放级别 15、水平位置 19144 和垂直位置 9524 的切片所覆盖的矩形。此处顶点的坐标为经纬度

POLYGON ((30.322265625 59.955010262062061, 30.322265625 59.949509172252277, 30.333251953125 59.949509172252277, 30.333251953125 59.955010262062061, 30.322265625 59.955010262062061))

标量 SQL 函数代码如下

图 1 切片边界几何图形创建列表

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

让我们依靠这种情况：如果在一个缩放级别上对象与切片没有交集，那么在下一个缩放级别上，与该切片匹配的更大编号的四个切片也没有交集。因此，当遍历切片和缩放级别时，只有当前一个缩放级别的切片与对象几何图形有交集时，才会进入下一个缩放级别进行交集检测。这有助于避免方法 2 中进行的过度检查。

• 方法 2

实际上，这是最坏的情况。因此，对于每个缩放级别，对于每个对象，通过使用函数 GEOMETRY::STEnvelope() 检测对象形状的边界框，并在该边界框内的所有切片中执行交集检测来收集切片子集。对于具有巨大面积或非常长长度的对象形状，这会严重影响性能。因为对于高缩放级别，它会导致检查大量切片。

• 方法 3

此方法应选择用于穿过大陆的复杂地理线。因此它适用于 LineString 几何类型。

对于边界框内瓦片边界网格中的每个对象，会形成一个单一值的几何图形集合。然后通过函数 GEOMETRY::STInersection 收集一组点。形状与网格几何图形的交点告诉我们要渲染哪些瓦片。因此，对于集合中的每个点，它会收集点之间的两个瓦片，并将瓦片位置编号添加到结果瓦片列表中。瓦片边界网格是一组垂直和水平线，它们穿过形状边界框中的瓦片边界。这比检查几何图形边界矩形区域中的每个瓦片更有效。

我们将公开第一种方法。

对于具有正方形的几何对象，例如 Polygon 和 Multipolygon，要检查交集的切片集被限制在覆盖对象形状的矩形区域的角切片中。

我们使用 CLR 函数 GEOMETRY::STEnvelope 获取对象形状边界区域的对角点，然后检测角切片位置。

对于几何图形形状类型为 Point 的对象，作为边界区域以检测与切片的交集，我们使用图标覆盖矩形区域。为了获得此图标覆盖矩形区域，编写了 CLR 函数 GetImageBound()。此函数在 GoogleProjection 类中。这里还有将经纬度转换为像素位置的方法。因为矩形区域的角点坐标表示为经纬度，我们需要将它们转换为切片位置编号。下一步，我们检测覆盖对象边界区域的切片子集，为此我们使用将坐标转换为指定缩放级别切片位置的函数。通过经纬度获取 X 和 Y 切片位置，可以使用已实现的 SQL CLR 函数或以下原生 SQL 函数

tile.GetXTilePos((@Longitude FLOAT, @Zoom INT)
tile.GetYTilePos((@Latitude FLOAT, @Zoom INT)

在我们获得角点切片位置后，通过函数 tile.fn_FetchGeometryTilesZoomDepth() 检查它们之间的所有切片与对象几何图形或图标覆盖区域的交集。

获取经度和缩放级别的 X 切片位置的 SQL 函数

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 切片位置的 SQL 函数

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 FLOAT
  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 中添加记录，并对当前切片下方的下一缩放级别的四个切片递归调用相同函数，并减小缩放深度。

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.Table，因为切片集将是唯一的。为了形成与不同对象有多个交集的切片，需要进行切片连接或重叠。因此，对于每个对象，将在分组过程中逐个创建切片并叠加在一起。

函数 tile.fn_GetTilesByTileNumZoomDepth 按照我们之前提到的，使用参数 @Depth 中指定的缩放深度执行交集检测。例如，如果参数 @Zoom = 1 和 @Depth = 1 的起始缩放级别为 1，则检查缩放级别 2 的切片，如果发生交集，则检查缩放级别 3 的四个切片。检查这些切片是必要的，因为它们覆盖了前一个缩放级别的切片。为了进行正确的交集检测，最好将 GEOMETRY 转换为 GEOGRAPHY。这很重要，因为当我们使用 GEOMETRY::STIntersects() 函数对 SRID = 4326 的 GEOGRAPHY 进行操作时，系统会知道坐标位于球体上，并在球体上进行交集计算。

折线对象形状示例

假设我们需要选择从圣彼得堡市中心到莫斯科市中心的折线切片。长度约为 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 个切片。下表显示了此形状按缩放级别分布的切片数量

图 2 显示了表示所构建折线的切片

图 2 - 切片 3.4.2 和 4.9.4

对于每个缩放级别的切片子集都会被收集，并且每个切片都会被检查是否存在交集。在第 4 和第 5 缩放级别，使用对象几何图形的函数 GEOMETRY::STEnvelope() 计算出的矩形区域内的切片数量并不大。请注意，在第 4 缩放级别只有 256 个切片。然而，在第 17 缩放级别，要检查的切片数量相当多。方法 1 有助于排除不必要的检查。对于几何图形为 POINT 的对象，每种方法都具有相同的效果。

交集检测

GEOMETRY::STIntersects() 函数有时在使用 GEOMETRY 类型而不是 GEOGRAPHY 时无法检测到交集。因为 GEOMETRY 函数使用 正交坐标 的公式，而 SRID = 4326 的 GEOGRAPHY 则使用球面坐标的公式。经纬度不是正交坐标。因此，为了进行精确的交集检测，正如我们前面提到的，我们需要将 GEOMETRY 值转换为 GEOGRAPHY。GEOGRAPHY 受多边形圆中点序列方向的限制。外部圆的坐标应按逆时针顺序排列。内部圆（孔）的坐标按顺时针顺序枚举。为避免转换错误，我们应使用函数 GEOMETRY::MakeValid() 和 GEOMETRY::STUnion() 将 GEOMETRY 转换为 GEOGRAPHY，然后获取多边形圆中正确的坐标序列。创建 GEOGRAPHY 时，将 SRID 设置为 4326——它告诉系统所有坐标都是球面坐标。(http://en.wikipedia.org/wiki/World_Geodetic_System)

目前，我们已将要渲染的切片列表存储在数据库表中，包含 Z、X、Y 列；渲染过程可以通过名为 mapnik 的渲染程序进行 (https://github.com/mapnik/mapnik/wiki/XMLConfigReference)。该程序的目的是从源地理空间信息生成切片。使用 mapnik 进行渲染配置需要一些工作，包括以下几个步骤

• 声明对象图层样式；
• 配置数据源；
• 设置切片列表源以避免为整个区域生成切片。

然而，在 SQL Server 引擎内部生成切片是可能的。为此，我们需要编写一些 SQL CLR 函数来操作几何数据类型。我们需要的主要函数是

• tile.IconTile(),
• tile.ShapeTile(),
• tile.TileAgg(),
• tile.SaveToFolderByZoomXY().

用于存储切片图像的表

图 2 显示了这些表。在这些表中，我们将存储要渲染的切片编号。字段 Data 的目的是以 PNG 格式存储切片图片。切片图片包含对象形状和图标的绘制。在数据库中操作大量图片会影响性能。对于此任务，更合适的情况是为每个对象从切片列表中在数据库外部形成切片，并按切片编号组合图片。最后，图片存储到文件系统。

图 3 – 存储切片图片的两张表

在特定位置的切片上绘制图标

让我们观察在特定缩放级别将图标定位到切片上的逻辑。这里我们讨论 POINT 几何类型。

我们有对象的经纬度和当前缩放级别要渲染的切片列表。切片列表的填充已在前面描述。解析图标在切片上的位置包括几个步骤

1. 将对象的经纬度转换为绝对像素编号；
2. 对于切片编号列表中的每个切片，对于当前缩放级别，计算切片左上角像素的像素编号。切片的左上角像素是通过将切片编号乘以 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>
/// Paste icon on 
/// </summary>
/// <param name="zoom"> zoom level number</param>
/// <param name="Lon">Longitude</param>
/// <param name="Lat">Latitude</param>
/// <param name="xTile">X tile number</param>
/// <param name="yTile">Y tile number</param>
/// <param name="iconImage">Icon bitmap</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 聚合 (TileAgg) 来实现瓦片分组。我们使用此聚合将不同对象但相同瓦片编号的瓦片图像（带有图标）组合成一个单一值。这里我们有一个缺点，对于每个对象，我们在数据库表中都有一个带有二进制字段的瓦片图像数据记录。这太多了。更有效的方法是只保留一个瓦片实例，并逐步绘制所有与它相交的对象图标。在这种情况下，内存消耗比我们使用瓦片分组时要少。尽管我们希望利用分组。

下一个过程将填充包含切片位置和绘制了对象图标的切片图像的表

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.Shape 表，其中存储了对象坐标、对象类型 ID 和图标图像的 ID，该图像存储在 tile.Image 表的 Binary 格式字段中。

下面的脚本是一个示例，将位置为 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)

图 4 – 由函数 tile.IconTile 形成的切片

在切片上绘制几何形状

对于折线几何形状 (POLYLINE, MULTIPOLYLINE)，我们获取瓦片边界与对象形状几何图形的几何交集，从而排除瓦片外部的部分。确定轮廓和填充区域的算法可应用于具有正方形的几何图形，而不适用于折线。因此，它适用于包含 POLYGON 或 MULTIPOLYGON 的 POLYGON, MULTIPOLYGON, GEOMETRYCOLLECTION。此算法的逻辑在 ShapeToTileRendering 类中实现，并包括以下步骤

1. 几何图形的坐标（经纬度）根据缩放级别通过将经纬度转换为像素坐标 PSG3857（Google 投影）的公式转换为像素位置编号，然后减去要渲染瓦片的左上角像素位置。假设在此步骤我们得到几何图形 (A)。此操作在函数 ConvertToZoomedPixelZeroedByTileGeometry(poly, Z, X, Y) 中实现；
2. 这里我们在零位置的切片像素坐标中创建几何图形，但根据缩放级别。这是几何图形 (B)；
3. 通过几何图形 (A) 和 (B) 的交集，我们得到了几何图形 (C)。这显然是 Polygon 或 Multipolygon。到目前为止，我们已经有了一个要填充的区域；
4. 形成几何图形 (D) 以移除几何图形 (C) 中穿过瓦片边界的线。因此，我们创建了几何图形 (C) 轮廓和瓦片几何图形 (B) 轮廓的交集，这意味着我们使用 Polyline 或 Multipolyline 形状类型而不是 Polygon 或 Multipolygon；
5. 几何图形 (E) 实际上是要绘制的轮廓。我们通过从 (C) 几何图形的轮廓中减去 (D) 来获得此轮廓。我们使用函数 GEOMETRY::STDifference 进行减法；
6. 我们在 (C) 几何图形内部区域的瓦片上进行填充——我们已经绘制了形状；
7. 我们在瓦片上绘制 (E) 几何图形 - 轮廓，并指定描边宽度，请注意在第 x 步已经排除了与瓦片边界的交集；
8. 对当前对象的所有瓦片重复步骤 1 到 7，并将它们保存到 Tile.TileOverlap 表中。

让我们通过示例更仔细地研究这些步骤。我们将为缩放级别 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 度，那么我们得到一个完整的圆，在其他情况下，它就像一个饼图或扇形。结果函数返回一个菱形的多边形几何图形

SELECT [tile].[fn_GetCircleSegment](30.3277587890625, 59.952259717159905,0,360,440,90)

要使用 Chrome 尝试此函数的各种参数值，请访问此实时示例。

SQL 函数代码如下

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

为了性能比较，我们实现了用于形成圆弧段几何图形的 CLR 函数和原生标量 SQL 函数。尽管它们在资源消耗方面没有明显差异。因此，这是 SQL CLR 函数的代码

/// <summary>
/// Builds tile geometry
/// </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).ToString()

结果如下

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

以上脚本的结果列在下表 2 的右侧两列中

现在我们准备好从获得的像素坐标创建切片轮廓和我们的菱形几何图形的交集几何图形

SELECT GEOMETRY::STGeomFromText('LINESTRING(4900936 2438400, 4901048 2438400, 4901120 2438328,
 4901120 2438216, 4901048 2438144, 4900936 2438144, 4900864 2438216, 4900864 2438328, 4900936 2438400
)',0)

在第 3 步形成的 (C) 几何图形是一个要填充的区域。它的轮廓在图 5 中以绿色绘制。

图 5 – 瓦片几何图形和对象几何图形的合并

我们对 (D) 几何图形并不是很感兴趣，尽管我们需要找出要在切片上绘制的作为我们的菱形轮廓的线条集合。这组线条，正如第 5 步所描述的，是通过从 (C) 几何图形的轮廓中减去 (D) 几何图形获得的。这组线条在图 5 中以蓝色绘制，可以用以下方式表示

SELECT GEOMETRY::STGeomFromText
 ('MULTILINESTRING(
  (4901048     2438400,      4901120       2438328),
  (4901120     2438216,      4901048       2438144),
  (4900936     2438144,      4900864       2438216),
  (4900864     2438328,      4900936       2438400))',
0)

图 6 – 要绘制的轮廓线集合（蓝线）

以下脚本创建 PNG 格式的切片图像，并将其保存到指定路径的文件系统。不存在的文件夹将根据切片存储结构创建

Z/X/Y.png

其中 Z – 缩放文件夹，X – 以列切片号命名的文件夹，Y.png – 以行切片号命名的文件。

DECLARE @bbox GEOMETRY
DECLARE @rhomb GEOMETRY
DECLARE @image VARBINARY(MAX)
SELECT @bbox = [tile].[GetTileBounds](15,19144,9524), 
       @rhomb = [tile].[fn_GetCircleSegment](30.3277587890625, 59.952259717159905,0,360,440,90)
SET @image = [tile].[ShapeTile]( @rhomb,15,19144,9524,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19144,9524)
SET @image = [tile].[ShapeTile]( @rhomb,15,19143,9524,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19143,9524)
SET @image = [tile].[ShapeTile]( @rhomb,15,19145,9524,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19145,9524)
SET @image = [tile].[ShapeTile]( @rhomb,15,19144,9523,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19144,9523)
SET @image = [tile].[ShapeTile]( @rhomb,15,19144,9525,'4400B050','9601B41E',3)
SELECT[tile].[SaveToFolderByZoomXY](@image,'d:/tiles',15,19144,9525)

通过上述脚本创建的 PNG 图像文件如图 7 所示。

图 7 - 菱形通过地理瓦片表示

/// <summary>
/// Proceed rendering of the geometry part on a tile
/// </summary>
/// <param name="shape">Full object geometry in degree coordinates (longitude, latitude)</param>
/// <param name="X">X tile number position</param>
/// <param name="Y">Y tile number position</param>
/// <param name="Zoom">number of zoom level</param>
/// <param name="argbFill">filling color in ARGB format</param>
/// <param name="argbStroke">color of a contour</param>
/// <param name="strokeWidth">contour width</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));
}

private void PasteShapeOnTile(Color fillcolor, Color strokecolor, int width, List<SqlGeometry> geom)
{
  SqlGeometry shape = geom[0];
  int geomnum = (int) shape.STNumGeometries();
  SqlGeometry stroke = null;
  SqlGeometry ring;
  int intnum;
  if (geom != null)
  switch (GetOpenGisGeometryType(shape))
  {
    case OpenGisGeometryType.LineString:
    case OpenGisGeometryType.MultiLineString:
         DrawMultiLineStringBordered2(shape, fillcolor, strokecolor, width, 1);
         break;
    case OpenGisGeometryType.Polygon:
         intnum = (int) shape.STNumInteriorRing();
         ring = shape.STExteriorRing();
         // 1. Draw polygon without inner holes
         FillPolygonOnTile(fillcolor, ring.ToPointsArray());
         // 2. Draw inner circles
         if (geomnum >= 1) stroke = geom[1];
         for (int i = 1; i <= intnum; i++)
         {
           FillTransparentPolygonOnTile(shape.STInteriorRingN(i).ToPointsArray());
         }
         // 3. Draw contour
         if (geom.Count > 1)
         {
           stroke = geom[1];
           DrawContourOnTile(stroke, strokecolor, width);
         }
         break;
    case OpenGisGeometryType.MultiPolygon:
         break;
  }
}

/// <summary>
/// Cuts shape by tile frame
/// </summary>
/// <param name="poly"></param>
/// <param name="X"></param>
/// <param name="Y"></param>
/// <param name="Z"></param>
/// <returns></returns>
private List<SqlGeometry> CutPolygonByZoomedPixelZeroTile(SqlGeometry poly, int X, int Y, int Z)
{
  return CutZoomedPixelPolygonByZeroTile(_parser.ConvertToZoomedPixelZeroedByTileGeometry(poly,Z,X,Y);
}

/// <summary
///  Proceeds conversion from degree coordinates to pixel and shifts them to ///zero tile       
/// </summary>
/// <param name="shape"></param>
/// <param name="zoom"></param>
/// <param name="tileX"></param>
/// <param name="tileY"></param>
/// <returns></returns>
public SqlGeometry ConvertToZoomedPixelZeroedByTileGeometry(SqlGeometry shape,int zoom,
         int tileX,int tileY)
{
  return CreateGeometryFromZoomedPixelInfo
          (ConvertToGeometryZoomedPixelsZeroTileShiftedInfo(
            GetGeometryInfo(shape), zoom, tileX, tileY));
}

/// <summary>
/// Converts from degree coordinates to pixels
/// Shifts to tile with position numbers equals to zero
/// </summary>
/// <param name="info"></param>
/// <param name="zoom"></param>
/// <param name="x"></param>
/// <param name="y"></param>
/// <returns>Geometry with pixel position numbers as coordinates</returns>
private GeometryZoomedPixelsInfo ConvertToGeometryZoomedPixelsZeroTileShiftedInfo
(GeometryInstanceInfo info, int zoom, int x, int y)
{
  int tilezeroshiftX = x*TILE_SIZE;
  int tilezeroshiftY = y*TILE_SIZE;
  var result = new GeometryZoomedPixelsInfo();
  var pixelCoordsListList = new List<List<GeometryPixelCoords>>();
  var geomPixCoordsList = new List<GeometryPixelCoords>();
  var coords = new GeometryPixelCoords {InnerRing = false};
  OpenGisGeometryType type = info.ShapeType;
  result.ShapeType = type;
  switch (type)
  {
    case OpenGisGeometryType.Point:
         PointF[] geopoints = info.Points[0][0].PointList;
         coords.PixelCoordList = new[]
           {new Point{X = _conv.FromLongitudeToXPixel(geopoints[0].X, zoom) - tilezeroshiftX,
                      Y = _conv.FromLatitudeToYPixel(geopoints[0].Y, zoom) - tilezeroshiftY }
           };
         geomPixCoordsList.Add(coords);
         pixelCoordsListList.Add(geomPixCoordsList);
         break;
    case OpenGisGeometryType.LineString:
         coords.PixelCoordList = 
            GetPixelCoordsShifted(
              info.Points[0][0].PointList,
              zoom,
              tilezeroshiftX,
              tilezeroshiftY);
         geomPixCoordsList.Add(coords);
         pixelCoordsListList.Add(geomPixCoordsList);
         break;
    case OpenGisGeometryType.Polygon:
         foreach (var list in info.Points)
           foreach (GeometryPointSequence pointseq in list)
           {
             coords.PixelCoordList = 
               GetPixelCoordsShifted(pointseq.PointList, zoom, tilezeroshiftX, tilezeroshiftY);
               coords.InnerRing = pointseq.InnerRing;
               geomPixCoordsList.Add(coords);
           }
         pixelCoordsListList.Add(geomPixCoordsList);
         break;
    case OpenGisGeometryType.MultiPoint:
    case OpenGisGeometryType.MultiLineString:
    case OpenGisGeometryType.MultiPolygon:
         pixelCoordsListList = 
         GetGeometryPixelCoordsShifted(info.Points, zoom, tilezeroshiftX, tilezeroshiftY);
         break;
    case OpenGisGeometryType.GeometryCollection:
         GeometryInstanceInfo[] geomColl = info.GeometryInstanceInfoCollection;
         int n = info.GeometryInstanceInfoCollection.Length;
         var geomPixZoomInfoCollection = new GeometryZoomedPixelsInfo[n];
         for (int i = 0; i < n; i++)
         {
             var geom = new GeometryZoomedPixelsInfo();
             geom.ShapeType = geomColl[i].ShapeType;
             geom.Points = 
             GetGeometryPixelCoordsShifted(geomColl[i].Points, zoom, tilezeroshiftX, tilezeroshiftY);
             geomPixZoomInfoCollection[i] = geom;
         }
         result.GeometryInstanceInfoCollection = geomPixZoomInfoCollection;
         break;
    }
    if (type != OpenGisGeometryType.GeometryCollection) result.Points = pixelCoordsListList;
  return result;
}

/// <summary>
/// Cuts already pixeled geometry by the frame of a tile
/// </summary>
/// <param name="poly"></param>
/// <param name="X"></param>
/// <param name="Y"></param>
/// <returns></returns>
private List<SqlGeometry> CutZoomedPixelPolygonByZeroTile(SqlGeometry poly, int X, int Y)
{
   List<SqlGeometry> result = new List<SqlGeometry>();
   SqlGeometry stroke = null;
   SqlGeometry contour;
   SqlGeometry tileLineString;
   SqlGeometry tobecut;
   SqlGeometry tile = _conv.GetTilePixelBound(0, 0, 1);
   var tiled = poly.STIntersection(tile);
   result.Add(tiled);
   switch (GetOpenGisGeometryType(tiled))
   {
      case OpenGisGeometryType.Polygon:
           // Create contour of polygon and inner rings as MULTILINESTRING
           contour = PolygonToMultiLineString(tiled);
           // remove cut lines of a geometry by tile frame
           tileLineString = tile.ToLineString();
           tobecut = contour.STIntersection(tileLineString);
           stroke = contour.STDifference(tobecut);
              break;
      case OpenGisGeometryType.MultiPolygon:
          // Create contour of multilinestring and inner rings as MULTILINESTRING
          contour = MultiPolygonToMultiLineString(tiled);
          // remove cutting lines by tile border
          tileLineString = tile.ToLineString();
          tobecut = contour.STIntersection(tileLineString);
          stroke = contour.STDifference(tobecut);
          break;
   }
   result.Add(stroke);
   return result;
}

[SqlFunction]
public static SqlBinary 
ShapeTile(SqlGeometry shape, SqlInt32 zoom,  SqlInt32 xTile, SqlInt32 yTile,
          SqlString argbFill,SqlString argbStroke,SqlInt32 strokeWidth)
{
  SqlBinary result = null;
  using (ShapeToTileRendering paster = new ShapeToTileRendering())
  {
    using (MemoryStream ms = new MemoryStream())
    {
      try
      {
        paster.DrawPartObjectShapeOnTile
          (shape,
           (int) xTile,
           (int) yTile,
           (int) zoom,
           argbFill.ToString(),
           rgbStroke.ToString(),
           (int) strokeWidth);
        result = paster.GetBytes();
      }
      catch (System.Exception ex)
      {
        string innerMessage = ex.InnerException.Message;
        throw 
          new Exception(
            string.Format("zoom: {1}; X:{2}; Y:{3} {0} , inner: {4}",
                          shape, zoom, xTile,yTile, innerMessage));
      }
      return result;
    }
  }
}

CREATE ASSEMBLY [Microsoft.SqlServer.Types]
AUTHORIZATION [dbo]
FROM 'd:\SQLCLR\BIN\TileRendering\Microsoft.SqlServer.Types.dll'
WITH PERMISSION_SET = UNSAFE
GO
CREATE ASSEMBLY [System.Drawing]
AUTHORIZATION [dbo]
FROM 'd:\SQLCLR\BIN\TileRendering\ System.Drawing.dll'
WITH PERMISSION_SET = UNSAFE
GO
CREATE ASSEMBLY [TileRendering]
AUTHORIZATION [dbo]
FROM 'd:\SQLCLR\BIN\TileRendering\TileRendering.dll'
WITH PERMISSION_SET = UNSAFE
GO
CREATE ASSEMBLY nQuant.Core
FROM 'd:\SQLCLR\BIN\TileRendering\ nQuant.Core.dll'
WITH PERMISSION_SET = UNSAFE
GO
CREATE ASSEMBLY SqlBitmapOperation
FROM 'd:\SQLCLR\BIN\TileRendering\SqlBitmapOperation.dll'
WITH PERMISSION_SET = UNSAFE
GO

ALTER DATABASE [dataBaseName] SET TRUSTWORTHY ON;

CREATE AGGREGATE [tile].[TileAgg]
(@Value [varbinary](max))
RETURNS[varbinary](max)
EXTERNAL NAME [SqlBitmapOperation].[TileAgg]
GO
CREATE AGGREGATE [tile].[IconTileAgg]
(@Value [varbinary](max), @PixelX [int], @PixelY [int])
RETURNS[varbinary](max)
EXTERNAL NAME [SqlBitmapOperation].[IconTileAgg]
GO
CREATE FUNCTION [tile].[IconTile]
(@image [varbinary](max), @zoom [int], @Lon [float], @Lat [float], @xTile [int], @yTile [int], @scale [float])
RETURNS [varbinary](max) WITH EXECUTE AS CALLER
AS 
EXTERNAL NAME [SqlBitmapOperation].[BitmapFunctions].[IconTile]
GO
--ShapeTile(SqlGeometry shape, SqlInt32 zoom,  SqlInt32 xTile, SqlInt32 yTile, SqlString argbFill,SqlString argbStroke,SqlInt32 strokeWidth)
CREATE FUNCTION [tile].[ShapeTile](@shape GEOMETRY, @zoom [int], @xTile [int], @yTile [int], @argbFill NVARCHAR(10),@argbStroke NVARCHAR(10), @strokeWidth INT)
RETURNS [varbinary](max) WITH EXECUTE AS CALLER
AS 
EXTERNAL NAME [SqlBitmapOperation].[BitmapFunctions].[ShapeTile]
GO
--SaveToFolderByZoomXY(SqlBinary image, SqlString rootFolderPath, SqlInt32 Zoom, SqlInt32 X,SqlInt32 Y)
CREATE FUNCTION tile.SaveToFolderByZoomXY(@image VARBINARY(MAX),@rootFolderPat NVARCHAR(512) , @Zoom [int], @xTile [int], @yTile [int])
RETURNS BIT WITH EXECUTE AS CALLER
AS 
EXTERNAL NAME [SqlBitmapOperation].[BitmapFunctions].[SaveToFolderByZoomXY]
GO

void PasteShapeOnTile(Color fillcolor,Color strokecolor, int width, List<SqlGeometry> geom)
{ 
  SqlGeometry shape = geom[0];
  int geomnum = (int)shape.STNumGeometries();
  SqlGeometry stroke = null;           
  SqlGeometry ring;
  int intnum;
  if  (geom != null)
  switch (GetOpenGisGeometryType(shape))
  {
    case OpenGisGeometryType.LineString:
    case OpenGisGeometryType.MultiLineString:
         DrawMultiLineStringBordered2(shape, fillcolor, strokecolor, width, 1);
         break;
    case OpenGisGeometryType.Polygon:                
         intnum = (int)shape.STNumInteriorRing();                   
         ring = shape.STExteriorRing();
         // 1. Draw polygon without inner rings
         FillPolygonOnTile(fillcolor, ring.ToPointsArray());
         // 2. Draw inner rings
         if (geomnum >= 1) stroke = geom[1];  
         for (int i = 1; i <= intnum; i++)
         {
           FillTransparentPolygonOnTile(shape.STInteriorRingN(i).ToPointsArray());
         }
         // 3. Draw contour
         if (geom.Count > 1)
         {
           stroke = geom[1];
           DrawContourOnTile(stroke, strokecolor, width);
         }
         break;
    case OpenGisGeometryType.MultiPolygon: break;
  }
}

CREATE PROC tile.FillShapeTiles  
            @GeoData GEOMETRY, @fillArgb VARCHAR(20),@strokeArgb VARCHAR(20),
            @FolderPath NVARCHAR(20), 
            @EndZoom INT = 17, @StartZoom INT = 4, @Thickness INT = 2
AS BEGIN
  IF @EndZoom < @StartZoom OR @GeoData IS NULL RETURN
  INSERT INTO tile.tile (Zoom, X,Y,Data)
  SELECT t.Zoom,
         t.TileX AS X,
         t.TileY AS Y,
         tile.ShapeTile(@GeoData, t.Zoom, t.TileX, t.TileY, @fillArgb, @strokeArgb ,@Thickness) AS Data   FROM 
   (SELECT * FROM tile.fn_FetchGeometryTilesZoomDepth(@GeoData,@StartZoom, @EndZoom - @StartZoom))    
  SELECT tile.SaveToFolderByZoomXY(Data, @FolderPath, Zoom, X, Y) 
  FROM tile.Tile         
END

//------------------------------------------------------------------------------
// <copyright file="CSSqlAggregate.cs" company="Microsoft">
//     Copyright (c) Microsoft Corporation.  All rights reserved.
// </copyright>
//------------------------------------------------------------------------------
using System;
using System.Collections.Generic;
using System.Data.SqlClient;
using System.Data.SqlTypes;
using Microsoft.SqlServer.Server;
using TileRendering;
using System.IO;
using System.Drawing;
using System.Drawing.Imaging;
[Serializable]
[SqlUserDefinedAggregate(Format.UserDefined, IsInvariantToDuplicates = true, IsInvariantToNulls = true, IsInvariantToOrder = false, IsNullIfEmpty = false, MaxByteSize = -1)]
public struct TileAgg : IBinarySerialize
{
  Bitmap _bitmap;
  ImageFormat _format;
  Graphics _graphics;
  ImageCodecInfo _codecInfo;
  const int TILE_SIZE = 256; 
  Bitmap GetInitialTile()
  {
    Bitmap DrawArea = new Bitmap(TILE_SIZE, TILE_SIZE);
    using (Graphics xGraph = Graphics.FromImage(DrawArea))
    {
      xGraph.FillRectangle(Brushes.Transparent, 0, 0, TILE_SIZE, TILE_SIZE);
      _graphics = Graphics.FromImage(DrawArea);
      return DrawArea;
    }
  }
#region [Aggregate artifacts]
  public void Init()
  {
    _codecInfo = GetEncoderInfo("image/png");
    _bitmap = GetInitialTile();
    DetectFormat();
  }

  public void Accumulate(SqlBinary Value)
  {
    using (MemoryStream ms = new MemoryStream())
    {
      ms.Write(Value.Value, 0, Value.Length);
      ms.Seek(0, SeekOrigin.Begin);
      ms.Position = 0;          
      PasteFromStreamImageToTile( ms);
    }
  }

  public void Merge(TileAgg Group)
  {
    PasteGroup(Group.Terminate());
  }

  public SqlBinary Terminate()
  {
    return GetBytes();
  }
#endregion [Aggregate artifacts]
  void PasteFromStreamImageToTile( Stream stream)
  {
    using (Bitmap iconImage = new Bitmap(stream, false))
    {
      DetectFormat();
      int width = iconImage.Width;
      int height = iconImage.Height;
      var area = new Rectangle(0, 0, width, height);
      CopyRegionIntoImage(iconImage,area, area);
    }
  }

  void CopyRegionIntoImage(Bitmap srcBitmap, Rectangle srcRegion, Rectangle destRegion)
  {
    _graphics.DrawImage(srcBitmap, destRegion, srcRegion, GraphicsUnit.Pixel);
    srcBitmap.Dispose();
  }
  
  void PasteGroup(SqlBinary Value)
  {
    using (MemoryStream ms = new MemoryStream())
    {
      ms.Write(Value.Value, 0, Value.Length);
      ms.Seek(0, SeekOrigin.Begin);
      ms.Position = 0;
      PasteTile(ms);
    }
  }

  void PasteTile(Stream stream)
  {
    Rectangle bounds = new Rectangle(0, 0, TILE_SIZE, TILE_SIZE);
    CopyRegionIntoImage(new Bitmap(stream), bounds, bounds);
  }

  byte[] GetBytes()
  {
    return _bitmap.ToByteArray(ImageFormat.Png);
  }

#region [IBinarySerialize]
  public void Read(BinaryReader reader)
  {
    _bitmap = new Bitmap(new MemoryStream(reader.ReadBytes((int)reader.BaseStream.Length)));
    DetectFormat();
  }

  public void Write(BinaryWriter writer)
  {
    EncoderParameters encodeParams = new EncoderParameters(1);
    encodeParams.Param[0] = new EncoderParameter(System.Drawing.Imaging.Encoder.Quality, 100);
    _bitmap.Save(writer.BaseStream, _codecInfo, encodeParams);
  }

#endregion [IBinarySerialize]
  /// <summary>
  /// Detects image format - very important
  /// </summary>
  void DetectFormat()
  {
    _format = _bitmap.GetImageFormat();
  }

  ImageCodecInfo GetEncoderInfo(string mimeType)
  {       
    string lookupKey = mimeType.ToLower();     
    ImageCodecInfo foundCodec = null;
    Dictionary<string, ImageCodecInfo> encoders = Encoders();
    if (encoders.ContainsKey(lookupKey))
      {
        foundCodec = encoders[lookupKey];
      }
    return foundCodec;
  }

  private Dictionary<string, ImageCodecInfo> Encoders()
  {
    Dictionary<string, ImageCodecInfo> encoders = new Dictionary<string, ImageCodecInfo>();  
    foreach (ImageCodecInfo codec in ImageCodecInfo.GetImageEncoders())
    {
      encoders.Add(codec.MimeType.ToLower(), codec);
    }
    return encoders;
  }
}

CREATE PROC tile.FillShapeTilesIntersection( @StartZoom INT, @EndZoom INT)
AS
BEGIN
DECLARE @Shape GEOMETRY   
DECLARE @CurrentZoom INT
DECLARE @ObjectTypeID INT
DECLARE @fillArgb NVARCHAR(10), @strokeArgb NVARCHAR(10)            
IF @ObjectTypeID IS NOT NULL
BEGIN        
  SET @CurrentZoom = @StartZoom                                        
  DECLARE shape_cursor CURSOR FOR 
  SELECT o.Shape, fillARGB, strokeARGB
  FROM tile.Shape o                                                   
  OPEN shape_cursor  
  FETCH NEXT FROM shape_cursor INTO @Shape, @fillArgb, @strokeArgb                                        WHILE @@FETCH_STATUS = 0
  BEGIN
    SET @CurrentZoom = @StartZoom
    WHILE @CurrentZoom  <= @EndZoom
    BEGIN
      INSERT INTO tile.tileOverlap (Zoom, X,Y,Data)
      SELECT t.Zoom, t.TileX AS X,t.TileY AS Y, 
             tile.ShapeTile(@Shape, t.Zoom, t.TileX, t.TileY, @fillArgb, @strokeArgb ,2) AS Data
      FROM (SELECT * FROM  tile.fn_FetchGeometryTiles(@Shape,@CurrentZoom)) t                                 SET @CurrentZoom = @CurrentZoom + 1
    END                              
    FETCH NEXT FROM shape_cursor INTO @Shape, @fillArgb, @strokeArgb 
  END
  CLOSE shape_cursor;
  DEALLOCATE shape_cursor;                     
  DELETE tile.TileOverlap
  END
END

INSERT INTO tile.Tile (Zoom, X, Y, Data)
SELECT Zoom, X, Y, tile.TileAgg(Data)
FROM tile.tileOverlap
GROUP BY Z, X, Y