STMorton: Creating a Morton Number Space Key Value for Grid Cells

Introduction

NOTE: This article is a copy of the original first written in 2010 from my old website.

Space curves are imaginary lines that systematically cover all tiles in an indexed space. The most well known in GIS is the Morton key. The utility of such a number can be seen though its implementing a form a space autocorrelation such that all geometries in an area are given a Morton number such that, when sorted, they appear next to each other in the sort order. Written to disk, such sorted geometries lie next to each other.

Function STMorton

 USE [GISDB] -- Change this to your database
 GO
 SET ANSI_NULLS ON
 SET QUOTED_IDENTIFIER ON
 GO
 -- =============================================
 -- Author:    <Author,Simon Greener>
 -- Create date: <Create Date, 7th August 2010>
 -- Description: <Description, 7th August 2010>
 -- =============================================
 CREATE FUNCTION [dbo].[STMorton] (@p_col INT, @p_row INT)
 RETURNS INT
 WITH EXECUTE AS CALLER,
      RETURNS NULL ON NULL INPUT
 AS
 BEGIN
   /* This procedure calculates the Morton number of a cell
      at the given row and column
      Written:  D.M. Mark, Jan 1984;
      Converted to Vax/VMS: July 1985
      Converted to Java: May 2007 (Simon Greener)
      Converted to PL/SQL: September 2009 (Simon Greener)
      Converted to TSQL: August 2010 (Simon Greener)
   */
 BEGIN
   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;
 END
 GO

Examples

Now for a simple test. (See dbo.MBR2Geometry )

 USE [GISDB] -- Change this to your database
 GO
 SELECT foo.MortonKey, foo.geom
   FROM (SELECT [GISDB].[dbo].STMorton( a.gcol, b.grow ) AS MortonKey,
          [GISDB].[dbo].MBR2GEOMETRY(a.gcol,b.grow,10.0,10.0,0) AS geom
           FROM (SELECT 0 + g.IntValue AS gcol FROM generate_series(0,7,1) AS g) AS a
                 CROSS APPLY
                (SELECT 0 + g.IntValue AS grow FROM generate_series(0,7,1) AS g) AS b
        ) foo
  ORDER BY mortonkey;

This is what it looks like:

gcolgrowmortonkeywktgeom
integerintegerintegertext
000POLYGON ((0 0,0 0,0 0,0 0,0 0))
101POLYGON ((1 0,1 0,1 0,1 0,1 0))
012POLYGON ((0 1,0 1,0 1,0 1,0 1))
113POLYGON ((1 1,1 1,1 1,1 1,1 1))
204POLYGON ((2 0,2 0,2 0,2 0,2 0))
305POLYGON ((3 0,3 0,3 0,3 0,3 0))
216POLYGON ((2 1,2 1,2 1,2 1,2 1))
317POLYGON ((3 1,3 1,3 1,3 1,3 1))
028POLYGON ((0 2,0 2,0 2,0 2,0 2))
129POLYGON ((1 2,1 2,1 2,1 2,1 2))
0310POLYGON ((0 3,0 3,0 3,0 3,0 3))
1311POLYGON ((1 3,1 3,1 3,1 3,1 3))
2212POLYGON ((2 2,2 2,2 2,2 2,2 2))
3213POLYGON ((3 2,3 2,3 2,3 2,3 2))
2314POLYGON ((2 3,2 3,2 3,2 3,2 3))
3315POLYGON ((3 3,3 3,3 3,3 3,3 3))
4016POLYGON ((4 0,4 0,4 0,4 0,4 0))
5017POLYGON ((5 0,5 0,5 0,5 0,5 0))
4118POLYGON ((4 1,4 1,4 1,4 1,4 1))

We can show the actual “curve” as a directed line through the morton numbers as follows. (The required function DumpPoints is described by following the link.)

 WITH morton_grid AS (
 SELECT foo.MortonKey, foo.geom
   FROM (SELECT [GISDB].[dbo].STMorton( a.gcol, b.grow ) AS MortonKey,
          [GISDB].[dbo].MBR2GEOMETRY(a.gcol,b.grow,10.0,10.0,0) AS geom
           FROM (SELECT 0 + g.IntValue AS gcol FROM [GISDB].[dbo].generate_series(0,7,1) AS g) AS a
                 CROSS APPLY
                (SELECT 0 + g.IntValue AS grow FROM [GISDB].[dbo].generate_series(0,7,1) AS g) AS b
        ) foo
 )
 SELECT geometry::STGeomFromText('LINESTRING(' +
        STUFF((SELECT ',' + a.coord
                  FROM (SELECT m.MortonKey, STR(e.[x],5,1) + ' ' +  STR(e.[y],5,1) AS coord
                  FROM morton_grid m
                 CROSS apply
                 [GISDB].[dbo].DumpPoints(m.geom.STCentroid()) e
                        ) a
                  ORDER BY ',' + STR(a.mortonkey,12,0)
               FOR XML PATH(''), TYPE, ROOT).VALUE('root[1]','nvarchar(max)'),1,1,'')
      + ')',0) AS geom;

Superimposing the curve over the grid of Morton keys reveals:

Conclusion

I hope this is of use to someone out there.