Space Filling Curve Functions For TSQL

Introduction

Warning: This article is long mainly because it includes the source code of the Hilbert, Morton and Peano key/curves.

I would assert that space filling curves have always been overlooked in their use within GIS.

Their existence was/is rarely taught at Universities, though at least in my day. Perhaps today it is different (but does ESRI support them?).

What they are can be seen from this Wikipedia Article and, also here.

Space Key values are generally created from the column/row references that index a grid and its cells. There are many Space Keys (thanks to Wikipedia); here are three:

• Hilbert Curve — First described by the German mathematician David Hilbert in 1891, as a variant of the space-filling Peano curves
• Peano Curve — Similar to Morton by Italian Mathematician Giuseppe Peano (1858 – 1932)
• Morton Curve— Named after Guy Macdonald Morton

These three implementations can be applied to any grid: it does not have to be part of Bing Maps’ tiling system.

Having mentioned Bing Maps, the other proprietary system is Google Maps.

What is the difference between the three I have mentioned?

Difference

This is what Wikipedia says about Morton Order/Code(Key):

“Morton space filling curve […] map[s] multidimensional data to one dimension while preserving locality of the data points.”

Which looks like this:

Morton Z-Order

The Peano and Hilbert curves have similar properties to, but are different from, Morton. The Peano curve is continuous curve that passes through every point of the unit square; it is described in this Wikipedia article. And the Hilbert Curve, a continuous fractal space-filling curve, is derived from the Peano, is described here.

Now, I’ll be frank, the theory and maths behind each of these is somewhat beyond me. But they have practical value in enabling us to take 2D data, index it using a space key, and then store the data in spatial order.

Tobler’s First Law of Geography is apposite here.

It is great to see PostGIS taking the lead with providing an embedded space curve capability when applying a SQL ORDER BY clause to a geometry/geography. I have done this for years manually and can only hope that what the PostGIS team has done will spread to other database.

Functions and Examples

I will present the three functions and after each function I will show its space curve over a common 16 col x 16 row grid.

But first we need to create the grid!

Base Grid

The following will build a 12 cell x 12 cell grid.

drop table if exists dbo.grid_cells;
go

with data as (
select 0.0 as rAngle,
       'LL' as anchor_point,
       10.0 as tileX,
       10.0 as tileY,
       32 as numTiles,
       geometry::Point(0,0,0) as origin -- LL
)
select 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;
go

This is what it looks like.

Hilbert Curve

Here is a translation of the C code on the Wikipedia article.

It only works for grids of size 4, 8, 16, 32, 64, etc. This is a geometric sequence that actually starts at 2 and is generated as follow:

2^1=2 — Cannot be used with Hilbert (or any space curve function)
2^2=4
2^3=8
2^4=16
2^5=32
2^6=64
2^7=128
2^n=xxx

CREATE FUNCTION [dbo].[STHilbert] 
(
  @p_n integer,  /* n is number of grid cells in a side -- must be same both sides*/
  @p_x integer,  /* Between 0 and @p_n */
  @p_y integer   /* Between 0 and @p_n */
)
RETURNS integer 
AS
/****f* SORT/STHilbert (2008)
 *  NAME
 *    STHilbert -- Function which creates a Hilbert (Space) Key from the supplied row and column reference.
 *  SYNOPSIS
 *    Function [dbo].[STHilbert] ( 
 *               @p_col bigint,
 *               @p_row int  
 *             )
 *     Returns int
 *  DESCRIPTION
 *    Function that creates a Hilbert Key from a row/col (grid) reference. 
 *    The generated value can be used to order/sort geometry objects.
 *  INPUTS
 *    @p_n   (int) -- The number of grid cells in a side ie 8 i.e. 8x8
 *    @p_col (int) -- Grid Column Reference beween 0 and @p_n
 *    @p_row (int) -- Grid Row Reference between 0 and @p_n.
 *  RESULT
 *    Hilbert_key  (int) - single integer Hilbert key.
 *  NOTES
 *    Works only for @p_n of 4, 8, 16, 32, 64 etc
 *  EXAMPLE
 *    SELECT dbo.STHilbert (16, 8, 8) as hKey;
 *    GO
 *
 *    hKey
 *    128
 *
 *    -- This example generates the hilbert key of a grid of cells.
 *    With data as (
 *      select x.IntValue as x, 
 *             y.IntValue as y
 *        from dbo.generate_series(0,3,1) as x,
 *             dbo.generate_series(0,3,1) as y
 *    )
 *    select geometry::STGeomFromText(
 *            'LINESTRING(' + STRING_AGG(dbo.STPointGeomAsText(point,12,12,12),',') WITHIN GROUP (order by hkey)+ ')',
 *            0).STAsText() as line
 *      from (select geometry::Point(d.x+0.5,d.y+0.5,0) as point, 
 *                   dbo.STHilbert(4,d.x,d.y) as hKey
 *             from data as d 
 *            ) as f
 *    go
 *
 *    LINESTRING (0.5 0.5, 1.5 0.5, 1.5 1.5, 0.5 1.5, 0.5 2.5, 0.5 3.5, 1.5 3.5, 1.5 2.5, 
 *                2.5 2.5, 2.5 3.5, 3.5 3.5, 3.5 2.5, 3.5 1.5, 2.5 1.5, 2.5 0.5, 3.5 0.5)
 *  AUTHOR
 *    Simon Greener
 *  HISTORY
 *    https://en.wikipedia.org/wiki/Hilbert_curve
 *    Simon Greener           - November 2020 - Original Coding for SQL Server.
 *  COPYRIGHT
 *    (c) 2008-2020 by TheSpatialDBAdvisor/Simon Greener
******/
BEGIN
  DECLARE
    @rx bigint,
    @ry bigint,
    @x  bigint,
    @y  bigint,
    @s  bigint,
    @t  bigint,
    @d  bigint;

  IF (@p_n is null or @p_x is null or @p_y is null)
    return null;

  SET @d  = 0;
  SET @x  = @p_x;
  SET @y  = @p_y;
  SET @s  = @p_n/2;
  While (@s>0)
  Begin
    SET @rx = case when (@x & @s) > 0 then 1 else 0 end; 
    SET @ry = case when (@y & @s) > 0 then 1 else 0 end; 
	SET @t = (3 * @rx) ^ @ry;
    SET @d  = @d + (@s * @s) * @t ;
    -- void rot(int n, int *x, int *y, int rx, int ry)
    -- Inline rot code 
    -- >>>>>>>
    If @ry = 0
    begin
      If @rx = 1 
      begin
        SET @X = @s - 1 - @X;
        SET @Y = @s - 1 - @Y;
      end;
      SET @T = @X;
      SET @X = @Y;
      SET @Y = @T;
    end;
    -- <<<<<<<
    SET @s = @s/2;
  END; /* while */
  Return @d;
END;
GO

All the functions need access to bitwise operators. Additionally access to the left and right shift bitwise operators is also needed. A separate TSQL implementation for these are available, but not shown here to limit the article’s length.

Now that the function has been defined, we can execute it, using praameters from the above common grid, as follows.

With data as (
  select x.IntValue as x, 
         y.IntValue as y
    from dbo.generate_series(0,15,1) as x,
         dbo.generate_series(0,15,1) as y
)
select geometry::STGeomFromText(
        'LINESTRING(' + STRING_AGG(dbo.STPointGeomAsText(point,12,12,12),',') WITHIN GROUP (order by hkey)+ ')',
        0).STBuffer(0.05) as line
  from (select geometry::Point(d.x+0.5,d.y+0.5,0) as point, 
               dbo.STHilbert(16,d.x,d.y) as hKey
         from data as d 
        ) as f
union all
select geom from dbo.grid_cells;
go

Visually this looks like:

The addition to the Point of 0.5 is to align, visually, with the grid. The STBuffer(0.05) is purely for visualisation.

Peano

The Peano key and curve works without the strict cell size requirement of the above STHilbert function. It will work regardless as the number of cells per side (though the resulting image may not be what is expected if the grid is not square).

The STPeano function, translated from some C source code provided to me by Dr David Mark, is as follows.

CREATE FUNCTION [dbo].[STPeano](
  @col bigint,
  @row bigint 
)
RETURNS bigint
As
/****f* SORT/STPeano (2008)
 *  NAME
 *    STPeano -- Function which creates a Peano (Space) Key from the supplied row and column reference.
 *  SYNOPSIS
 *    Function STPeano ( 
 *               @p_col int,
 *               @p_row int  
 *             )
 *     Returns int
 *  USAGE
 *  DESCRIPTION
 *    Function that creates a Peano Key from a row/col (grid) reference. 
 *    The generated value can be used to order/sort geometry objects.
 *  INPUTS
 *    @p_col      (int) - Grid Column Reference.
 *    @p_row      (int) - Grid Row Reference.
 *  RESULT
 *    Peano_key  (int) - single integer Peano key.
 *  EXAMPLE
 *    SELECT dbo.STPeano (10, 10) as pKey;
 *    GO
 *
 *    pKey
 *    204
 *
 *    SELECT [dbo].[STPeano](row.IntValue,col.IntValue) as pkey
 *      FROM [dbo].[generate_series](0,3,1) as col,
 *           [dbo].[generate_series](0,3,1) as row
 *     ORDER BY peano_key;
 *    GO
 *    
 *    peano_key
 *    0
 *    1
 *    2
 *    3
 *    4
 *    5
 *    6
 *    7
 *    8
 *    9
 *    10
 *    11
 *    12
 *    13
 *    14
 *    15
 *  AUTHOR
 *    Simon Greener
 *  HISTORY
 *    Professor David M. Mark - January  1984 - C;
 *    Simon Greener           - November 2020 - Original Coding for SQL Server.
 *  COPYRIGHT
 *    (c) 2008-2020 by TheSpatialDBAdvisor/Simon Greener
******/
Begin
  DECLARE
    @curveaddr bigint,
    @curvemask bigint,
    @i         int;

    /* Compute space curve address for point */
    /* PEANO: bit interleave col and row */

    SET @i = 16; -- only use low 16 bits of col and row 
    SET @curveaddr = 0; 
    SET @curvemask = 1;
    while (@i > 0)
    BEGIN
        set @curveaddr = @curveaddr | case when (@row & 1) <> 0 then @curvemask else 0 end;
        set @curvemask = [dbo].[shiftLeft](@curvemask,1); 
        set @row       = [dbo].[shiftRight](@row,1);
        set @curveaddr = @curveAddr | case when (@col & 1) <> 0 then @curvemask else 0 end;
        set @curvemask = [dbo].[shiftLeft](@curvemask,1);
        SET @col       = [dbo].[shiftRight](@col,1);
        set @i -= 1;
    END;
    RETURN @curveaddr;
END;
GO

Now that the function has been defined, we can execute it, using praameters from the above common grid, as follows.

With data as (
  select x.IntValue as x, 
         y.IntValue as y
    from dbo.generate_series(0,15,1) as x,
         dbo.generate_series(0,15,1) as y
)
select geometry::STGeomFromText(
        'LINESTRING(' + STRING_AGG(dbo.STPointGeomAsText(point,12,12,12),',') WITHIN GROUP (order by hkey)+ ')',
        0) as line
  from (select geometry::Point(d.x+0.5,d.y+0.5,0.5) as point, 
               dbo.STPeano(d.x,d.y) as hKey
         from data as d 
        ) as f
union all
select geom from dbo.grid_cells;
go

Visualising we get:

Morton

Finally, the Morton key which creates a space key and code from the row and column references.

The STMorton function, translated from some C source code provided to me by Dr David Mark, is as follows.

CREATE FUNCTION [dbo].[STMorton] 
(
  @p_col int, 
  @p_row int
)
Returns int
WITH EXECUTE AS CALLER, 
     RETURNS NULL ON NULL INPUT
AS
BEGIN
/****f* SORT/STMorton (2008)
 *  NAME
 *    STMorton -- Function which creates a Morton (Space) Key from the supplied row and column reference.
 *  SYNOPSIS
 *    Function STMorton ( 
 *               @p_col int,
 *               @p_row int  
 *             )
 *     Returns int
 *  DESCRIPTION
 *    Function that creates a Morton Key from a row/col (grid) reference. 
 *    The generated value can be used to order/sort geometry objects.
 *  INPUTS
 *    @p_col      (int) - Grid Column Reference.
 *    @p_row      (int) - Grid Row Reference.
 *  RESULT
 *    morton_key  (int) - single integer morton key.
 *  EXAMPLE
 *    SELECT STMorton (10, 10) as mKey;
 *     # mKey
 *     828
 *  AUTHOR
 *    Simon Greener
 *  HISTORY
 *    Professor David M. Mark - January  1984 - C;
 *    Simon Greener           - December 2011 - Original Coding for SQL Server.
 *  COPYRIGHT
 *    (c) 2008-2018 by TheSpatialDBAdvisor/Simon Greener
******/
  Declare 
    @row       int = abs(@p_row),
    @col       int = abs(@p_col),
    @key       int = 0,
    @level     int = 0,
    @left_bit  int,
    @right_bit int,
    @quadrant  int;
  Begin
    While ((@row>0) Or (@col>0))
    Begin
      /* Split off the row (left_bit) and column (right_bit) bits and then combine them to form a bit-pair representing the quadrant */
      Set @left_bit  = @row % 2;
      Set @right_bit = @col % 2;
      Set @quadrant  = @right_bit + 2*@left_bit;
      Set @key       = @key + round(@quadrant * power(2,2*@level), 0, 1);
      /* row, column, and level are then modified before the loop continues */
      If ( @row = 1 And @col = 1 )
      Begin 
        Set @row = 0; 
        Set @col = 0;
      End
      Else
      Begin
        Set @row = @row / 2;
        Set @col = @col / 2;
        Set @level = @level + 1;
      End;
    End;
  End;
  Return @key;
End;
GO

Now that the function has been defined, we can execute it, using parameters from the above common grid, as follows.

With data as (
  select x.IntValue as x, 
         y.IntValue as y
    from dbo.generate_series(0,15,1) as x,
         dbo.generate_series(0,15,1) as y
)
select geometry::STGeomFromText(
        'LINESTRING(' + STRING_AGG(dbo.STPointGeomAsText(point,12,12,12),',') WITHIN GROUP (order by hkey)+ ')',
        0).STBuffer(0.05) as line
  from (select geometry::Point(d.x+0.5,d.y+0.5,0.5) as point, 
               dbo.STMorton(d.x,d.y) as hKey
         from data as d 
        ) as f
union all
select geom from dbo.grid_cells;
go

Visually, the curve looks like the following.

Conclusion

The visualisation of the generation of the three space curves from their relevant functions shows how powerful they can be.

Tobler’s First Law of Geography is nobly implemented using these functions.

PostGIS is demonstrating this with its improvement of the ORDER BY in their product via use of Hilbert.

In the days before SSD, sorting and caching geometry objects on disk using a space curve can provide significant performance improvements.

Morton or Peano curves are best (practically) for the simple reason that you do not have to specific a strict number of grid cells, though that does not mean that Hilbert does not have its uses (cf PostGIS).

I’ve used Peano and Morton before over the years to achieve my needs at the time. They work, and I will continue to use them. I created a Hilbert function because of the PostGIS announcement and for this article.

This article can be edited if anyone finds a error, or can provide improved descriptions of space curve/keys and my functions.