Gridding a geometry or geography object (SQL Server Denali)

A common question that comes up in the various database forums (PostGIS, SQL Server, Oracle) is how to “grid” a linear or polygonal object.

By “grid” one means work out the square pixels (rectangular polygons) that cover or define a vector geometry.

Here is some SQL for doing this for SQL Server 2008 Denali. Denali has spatial aggregates built in which is why it was used for this blog article. It can be done for SQL Server 2008 R1 and R2 but it needs the Codeplex extensions to be installed and enabled.

The gridding method requires two functions Morton and REGULARGRIDXY. THe Morton key function is provided in another blog article.

The other required function, REGULARGRIDXY ensures that the gridding takes place according to a commonly defined MBR so that the grids generated across multiple geometries align. This function is provided here:

 CREATE FUNCTION [dbo].[REGULARGRIDXY]
 (
   @p_ll_x       FLOAT,
   @p_ll_y       FLOAT,
   @p_ur_x       FLOAT,
   @p_ur_y       FLOAT,
   @p_TileSize_X FLOAT,
   @p_TileSize_Y FLOAT,
   @p_srid       INT
 )
 RETURNS @TABLE TABLE
 (
   col   INT,
   ROW   INT,
   geom  geometry
 )
 AS
 BEGIN
    DECLARE
      @v_loCol INT,
      @v_hiCol INT,
      @v_loRow INT,
      @v_hiRow INT,
      @v_col   INT,
      @v_row   INT;
    BEGIN
      SET @v_loCol = FLOOR(   @p_LL_X / @p_TileSize_X );
      SET @v_hiCol = CEILING( @p_UR_X / @p_TileSize_X ) - 1;
      SET @v_loRow = FLOOR(   @p_LL_Y / @p_TileSize_Y );
      SET @v_hiRow = CEILING( @p_UR_Y / @p_TileSize_Y ) - 1;
      SET @v_col = @v_loCol;
      WHILE ( @v_col <= @v_hiCol )
      BEGIN
        SET @v_row = @v_loRow;
        WHILE ( @v_row <= @v_hiRow )
        BEGIN
          INSERT INTO @TABLE (col,ROW,geom)
          VALUES(@v_col, @v_row,
                 geometry::STGeomFromText('POLYGON((' +
                       LTRIM(STR((@v_col * @p_TileSize_X),24,12)) + ' ' + LTRIM(STR((@v_row * @p_TileSize_Y),24,12)) + ',' +
                       LTRIM(STR(((@v_col * @p_TileSize_X)+@p_TileSize_X),24,12)) + ' ' + LTRIM(STR((@v_row * @p_TileSize_Y),24,12)) + ',' +
                       LTRIM(STR(((@v_col * @p_TileSize_X)+@p_TileSize_X),24,12)) + ' ' + LTRIM(STR(((@v_row * @p_TileSize_Y)+@p_TileSize_Y),24,12)) + ',' +
                       LTRIM(STR((@v_col * @p_TileSize_X),24,12)) + ' ' + LTRIM(STR(((@v_row * @p_TileSize_Y)+@p_TileSize_Y),24,12)) + ',' +
                       LTRIM(STR((@v_col * @p_TileSize_X),24,12)) + ' ' + LTRIM(STR((@v_row * @p_TileSize_Y),24,12)) + '))',@p_srid)
                );
          SET @v_row = @v_row + 1;
        END;
        SET @v_col = @v_col + 1;
      END;
      RETURN;
    END;
 END
 GO

Here is the gridding method for a single geometry.

 -- The following should be changed to the database name in which you have installed the Morton and RegularGridXY functions
 --
 USE [GISDB]
 GO
 -- Single geometry processing
 --
 WITH geomQuery AS (
 SELECT g.geom.STEnvelope().STPointN(1).STX AS minx, g.geom.STEnvelope().STPointN(1).STY AS miny,
        g.geom.STEnvelope().STPointN(3).STX AS maxx, g.geom.STEnvelope().STPointN(3).STY AS maxy,
        g.geom, 0.050 AS gridX, 0.050 AS gridY, 0 AS loCol, 0 AS loRow
   FROM (SELECT a.geom.STBuffer(1.000).STSymDifference(a.geom.STBuffer(0.500)) AS geom
           FROM (SELECT geometry::STGeomFromText('MULTIPOINT((09.25 10.00),(10.75 10.00),(10.00 10.75),(10.00 9.25))',0) AS geom ) a
        ) g
 )
 SELECT f.mKey, f.col, f.ROW, f.geom
   FROM (SELECT [GISDB].[dbo].Morton((c.col - a.loCol),(c.ROW - a.loRow)) AS mKey,
                c.col, c.ROW,
          CASE WHEN UPPER(a.geom.STGeometryType()) IN ('POLYGON','MULTIPOLYGON')
               THEN a.geom.STIntersection(c.geom)
           ELSE a.geom
         END AS geom
           FROM geomQuery a
                CROSS APPLY
                [GISDB].[dbo].REGULARGRIDXY(a.minx, a.miny, a.maxx, a.maxy,a.gridX,a.gridY, a.geom.STSrid) c
          WHERE a.geom.STIntersects(c.geom) = 1
        ) f
  WHERE UPPER(f.geom.STGeometryType()) IN ('POLYGON','MULTIPOLYGON') /* Don't want point or line tiles */
  ORDER BY f.mKey;

The result of this looks like.

For multiple geometries, the following needs to be executed.

 -- Query for more than one geometry
 --
 WITH geomQuery AS (SELECT g.rid,
                           MIN(g.geom.STEnvelope().STPointN(1).STX) OVER (partition BY g.pid) AS minx,
                           MIN(g.geom.STEnvelope().STPointN(1).STY) OVER (partition BY g.pid) AS miny,
                           MAX(g.geom.STEnvelope().STPointN(3).STX) OVER (partition BY g.pid) AS maxx,
                           MAX(g.geom.STEnvelope().STPointN(3).STY) OVER (partition BY g.pid) AS maxy,
                           g.geom, 0.050 AS gridX, 0.050 AS gridY, 0 AS loCol, 0 AS loRow
                      FROM (SELECT 1 AS pid, a.rid, a.geom.STBuffer(1.000).STSymDifference(a.geom.STBuffer(0.750)) AS geom
                              FROM (SELECT 1 AS rid, geometry::STGeomFromText('POINT(09.50 10.00)',0) AS geom
                          UNION ALL SELECT 2 AS rid, geometry::STGeomFromText('POINT(10.50 10.00)',0) AS geom
              UNION ALL SELECT 3 AS rid, geometry::STGeomFromText('POINT(10.00 10.50)',0) AS geom
                          UNION ALL SELECT 4 AS rid, geometry::STGeomFromText('POINT(10.00 09.50)',0) AS geom ) a
              ) g
                   )
 SELECT ROW_NUMBER() OVER (ORDER BY f.col, f.ROW) AS tid,
        [GISDB].[dbo].Morton((f.col - f.loCol),(f.ROW - f.loRow)) AS mKey,
        f.col,
        f.ROW,
        COUNT(*) AS UnionedTileCount,
        geometry::UnionAggregate(f.geom) AS geom
   FROM (SELECT c.col, c.ROW, a.loCol, a.loRow,
            CASE WHEN UPPER(a.geom.STGeometryType()) IN ('POLYGON','MULTIPOLYGON')
               THEN a.geom.STIntersection(c.geom)
           ELSE a.geom
         END AS geom
           FROM geomQuery a
                CROSS APPLY
                [GISDB].[dbo].REGULARGRIDXY(a.minx,a.miny,a.maxx,a.maxy,a.gridX,a.gridY,a.geom.STSrid) c
          WHERE a.geom.STIntersects(c.geom) = 1
        ) f
  WHERE UPPER(f.geom.STGeometryType()) IN ('POLYGON','MULTIPOLYGON') /* Don't want point or line tiles */
  GROUP BY f.col, f.ROW, f.loCol, f.loRow
  ORDER BY 2;

That looks like this.

I hope this is of use to SQL Server Denali users