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;
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_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))'

I hope this is of use to someone.