Morton Key Function for PostgreSQL/PostGIS

In my article on the creation and use of a space keys for the spatial sorting of spatial data, I show how to implement a Morton key for PostgreSQL/PostGIS pl/pgSQL.

The function is as follows:

 CREATE OR REPLACE FUNCTION ST_Morton (p_col int8, p_row int8)
  RETURNS int8
 AS
 $$
 /*  this procedure calculates the Morton number of a cell
     at the given row and col[umn]
     Written:  D.M. Mark, Jan 1984;
     Converted to Vax/VMS: July 1985
     Converted to PostgreSQL, Simon Greener, 2010
 */
 DECLARE
    v_row          int8 := 0;
    v_col          int8 := 0;
    v_key          int8;
    v_level        int8;
    v_left_bit     int8;
    v_right_bit    int8;
    v_quadrant     int8;
 BEGIN
    v_row   := p_row;
    v_col   := p_col;
    v_key   := 0;
    v_level := 0;
    WHILE ((v_row>0) OR (v_col>0)) LOOP
      /* Split off the row (left_bit) and column (right_bit) bits and
      then combine them to form a bit-pair representing the quadrant */
      v_left_bit  := v_row % 2;
      v_right_bit := v_col % 2;
      v_quadrant  := v_right_bit + 2*v_left_bit;
      v_key       := v_key + ( v_quadrant << (2*v_level) );
      /* row, column, and level are then modified before the loop continues */
      v_row := v_row / 2;
      v_col := v_col / 2;
      v_level := v_level + 1;
    END LOOP;
    RETURN (v_key);
 END;
 $$
   LANGUAGE 'plpgsql' IMMUTABLE;

Now, let’s test it with some simple SQL:

 SELECT *
   FROM (SELECT a.gcol,
                b.grow,
                st_morton( a.gcol, b.grow ) AS MortonKey,
                ST_MakeBox2D(ST_MakePoint(a.gcol, b.grow)::geometry,
                             ST_MakePoint(a.gcol+100, b.grow+100)::geometry)::geometry AS geometry           FROM (SELECT 0 + g AS gcol FROM generate_series(0,7,1) AS g) AS a,
                (SELECT 0 + g AS grow FROM generate_series(0,7,1) AS g) AS b
         ) AS foo
  ORDER BY mortonkey, gcol;

Results

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

I hope this is of use to someone.