Geocoding Points Against a Grid of Cells/Tiles With/out Spatial Operators

Southern Tasmania

Let’s now construct an more realistic example that uses a rotated grid in a real place near home!

We will endeavor to discover the comparative performance of these two methods, though straight comparison depends on real situations in which each one solution may be best but the other not. (See conclusion.)

Now let’s use a different function to create an angled grid in a real world location.

/****f* TILING/STTileGeomByPoint (2008)
 *  NAME
 *    STTileGeomByPoint -- Creates mesh of tiles in geometric space anchored to a supplied origin.
 *  SYNOPSIS
 *    Function STTileGeomByPoint (
 *      @p_point    geometry,
 *      @p_origin   varchar(2),
 *      @p_numTileX integer,
 *      @p_numTileY integer,
 *      @p_TileX    float,
 *      @p_TileX    float,
 *      @p_rAngle   float,
 *      @p_AsPoint  bit
 *    )
 *    Returns @table table
 *    (
 *      col  Int,
 *      row  Int,
 *      geom geometry
 *    )
 *  DESCRIPTION
 *    This function generates a mesh (grid) of tiles anchored to the supplied anchor point.
 *    The mesh of tiles is controlled by five parameters:
 *      1. An anchor point
 *      2. Whether the anchor point is at the LowerLeft (LL), UpperRight (UR), UpperLeft (UL)L or LowerRight (LR) of the grid
 *      3. The number of tiles in X and Y direction;
 *      4. XY tile size in meters; 
 *      5. Optional rotation angle (around origin/anchor point)
 *  INPUTS
 *    @p_point  (geometry) -- Starting Point for grid (Upper Left)
 *    @p_origin   (varchar) -- Position of point wrt grid: LL,UL,LR,UR
 *    @p_numTileX (integer) -- Number of tiles in X (longitude) direction
 *    @p_numTileY (integer) -- Number of tiles in Y (latitude) direction
 *    @v_TileX      (float) -- Size of a Tile's X dimension in real world units along parallel of Latitude (ie X distance)
 *    @v_TileY      (float) -- Size of a Tile's Y dimension in real world units along meridian of Longitude (ie Y distance)
 *    @p_rAngle     (float) -- Optional rotation angle from North.
 *    @p_AsPoint      (bit) -- Return Tile as point or polygon
 *  RESULT
 *    A Table of the following is returned
 *    (
 *      col  Int      -- The column reference for a tile
 *      row  Int      -- The row reference for a tile
 *      geom geometry -- The polygon geometry covering the area of the Tile.
 *    )

Create and Store a Grid

So, let’s create, and store, an angled grid, with a Morton key as its primary key, and whose anchor point is the LL of the grid.

drop table if exists dbo.grid_cells 
go

with data as (
select -100.0 as rAngle,
       'LL'   as anchor_point,
       1000.0 as tileX,
       1000.0 as tileY,
       15     as numTiles,      -- 8x8 = 64
       64     as PointsPerTile, -- 8*8*64 = 4096
       geometry::Point(537038,5215350,28355) as origin -- LL
)
select [dbo].[STMorton](t.col,t.row) as morton_id, col, row, geom
  into dbo.grid_cells 
  from data as a
       cross apply
        [dbo].[STTileGeomByPoint] (
            a.origin,
            /*p_origin*/    a.anchor_point,
            /*@p_numTileX*/ a.numTiles,
            /*@p_numTileY*/ a.numTiles, /* 64 cells */
            /*@p_TileX   */ a.tileX,
            /*@p_TileY   */ a.tileY,
            /*@p_rAngle  */ a.rAngle,
            /*@p_AsPoint*/  0
         ) as t;

alter table dbo.grid_cells alter column morton_id int not null;
go

alter table dbo.grid_cells add constraint grid_cells_pk primary key clustered (morton_id) ;
go

CREATE SPATIAL INDEX [grid_cells_spx] ON [dbo].[grid_cells]([geom]) 
 USING  GEOMETRY_AUTO_GRID 
  WITH (BOUNDING_BOX =(200000, 5100000, 500000, 5500000))
    ON [PRIMARY]
GO

The generated grid looks like this.

Shows Grid Generated From TSQL Functions
Grid in Southern Tasmania.

Create And Store Sample Points – Fish, Rocks, Boats?

We can now will generate a set of points for our angled grid which can be used to compare the two methods.

Note that we cannot use the cells to directly create sample points because the grid is angled (it is easy if it is not). The approach the table creation code below uses is to take all the grid’s cells and merge them into a single object; then compute the single polygon’s envelope and then use it to generate a set of candidate points. Once the candidate points are created, merge them back to the merged single polygon keeping only those that fall within it.

drop table if exists dbo.grid_points;
go

with data as (
select -100.0 as rAngle,
       'LL'   as anchor_point,
       1000.0 as tileX,
       1000.0 as tileY,
       8      as numTiles,      -- 8x8 = 64
       96     as PointsPerTile, -- 8*8*96 = 6144
       geometry::Point(537038,5215350,28355) as origin -- LL
)
, mbr as (
-- Generate single tile to create area for creating points
select mbr.minx, mbr.maxx, mbr.miny, mbr.maxy
  from (select geometry::UnionAggregate(a.geom) as geom
          from dbo.grid_cells as a
       ) as t
       cross apply
       [dbo].[STGeometry2MBR] (t.geom) as mbr
) -- Now create points in this aggregate polygon (will remove that which are not needed later)
, points as (
SELECT gs.IntValue as point_id,
       geometry::Point(
          [dbo].[STRandomBetween](m.minx, m.maxx),
          [dbo].[STRandomBetween](m.miny, m.maxy),
          a.origin.STSrid
       ) as geom
  FROM data as a CROSS APPLY [dbo].[generate_series](1,a.numTiles * a.numTiles * a.PointsPerTile,1) as gs,
       mbr as m
)
SELECT f.point_id, f.geom
  INTO dbo.grid_points
  FROM (SELECT p.point_id, 
               row_number() over (partition by p.point_id order by p.point_id) as low_point_id,
               p.geom
          FROM points as p
               inner join
               dbo.grid_cells as t
               on (t.geom.STIntersects(p.geom) = 1) -- select only point that fall within single polygon
       ) as f
 WHERE low_point_id = 1
GO

alter table dbo.grid_points alter column point_id int not null;
go

alter table dbo.grid_points add constraint grid_points_pk primary key clustered (point_id) ;
go

CREATE SPATIAL INDEX [grid_points_spx] ON [dbo].[grid_points]([geom]) 
 USING  GEOMETRY_AUTO_GRID 
  WITH (BOUNDING_BOX =(200000, 5100000, 500000, 5500000))
    ON [PRIMARY]
GO

select count(*) from dbo.grid_points;
GO
-- 3263 out of 6144

These points can be visualised over the top of our grid as follow.

This image shows our generated points over the grid
Generated Points over Angled Grid

Note that there is a reduction in the number of points. They were generated for a square polygon over the angled grid; that mismatch causes the point reduction.

Comparison

We will now compare the two methods with our real world grid.

First the SQL Server Spatial (point in polygon) method.

Note that the result will include the two morton key values side by side, so we can visually compare the two morton keys (or we could filter as well).

-- Note, that the table has a spatial index.
select t.morton_id, 
       [dbo].[STMorton](t.col,t.row) as morton_point_id, 
       p.geom
  from dbo.grid_points as p
       inner join
       dbo.grid_cells as t
       on (t.geom.STContains(p.geom) = 1)
go

-- 3263 Rows in 3 seconds.

Next the direct reference approach.

with pointsWithRef as (
-- Create Morton Key from generated col/row references.
SELECT [dbo].[STMorton](r.col,r.row) as morton_point_id, r.col, r.row, f.geom
  FROM dbo.grid_points as f
       CROSS APPLY 
       [dbo].[STTileColRowRef4GeomPoint](
         f.geom,
         geometry::Point(537038,5215350,28355),
         1000.0,
         1000.0,
        -100.0
       ) as r
)
select t.morton_id, morton_point_id
  from pointsWithRef as p 
       inner join 
       dbo.grid_cells as t 
       on (t.morton_id = p.morton_point_id)
go

-- 3263 Rows in 7 seconds

Even with STContains’s best case use, SQL Server Spatial’s operator was still 2.3 times (7 sec / 3 secs) or 233% faster.

If we look at this from the point of view of speed / geometry we get:

  • 7 / 3263 = 0.002 seconds per feature, and
  • 3 / 3263 = 0.0009 seconds per feature.

Both are pretty quick.

But is the difference really that bad?