Hexagonal tiling for PostGIS

hexagon_tiling.sql

-- Hexagonal tiling for PostGIS
-------------------------------
-- Edward Abraham, Dragonfly Science
-- Use freely

-- Functions for generating tiles in a hexagonal tiling, from cartesian
-- coordinates. The tiling is made of hexagons defined by their 'width'
-- (the distance bewteen two parallel sides). The hexagon with index (0, 0)
-- is centered on the point x=0, y=0, and is oriented so that it points 
-- upwards. Hexagons in the same row have the same j-coordinate, and hexagons
-- along lines at 60 degrees anti-clockwise from the x-axis have the same
-- i-coordinate.

-- Example usage:

-- Find the tile (in a tiling of hexagons of width 15) that includes the 
-- point 200.0, 136.0
-- SELECT get_hexagon(200.0, 136.0, 15);


-- The algorithm is from the description given by James McNeill
-- http://playtechs.blogspot.co.nz/2007/04/hex-grids.html

-- I referred to the python implementation by Eric Gradman for
-- hexagon_i and hexagon_j
-- https://gist.github.com/583180 

-- Returns the i-index of the hexagon in a hexagonal tiling,
-- containing the cartesian point 'x', 'y'
-- Called with arguments i, j, width
CREATE OR REPLACE FUNCTION hexagon_i(FLOAT, FLOAT, FLOAT)
RETURNS INTEGER AS $$
    SELECT FLOOR((FLOOR($1/$3 - SQRT(3.0)*$2/$3) + FLOOR(2*$1/$3) + 2.0)/3.0)::INTEGER
$$ LANGUAGE SQL IMMUTABLE;

-- Returns the j-index of the hexagon in a hexagonal tiling,
-- containing the cartesian point 'x', 'y'
-- Called with arguments i, j, width
CREATE OR REPLACE FUNCTION hexagon_j(FLOAT, FLOAT, FLOAT)
RETURNS INTEGER AS $$
    SELECT FLOOR((FLOOR(($1 + SQRT(3.0)*$2)/$3) + FLOOR((-$1 + SQRT(3.0)*$2)/$3)+ 2.0)/3.0)::INTEGER
$$ LANGUAGE SQL IMMUTABLE;

-- Returns the hexagonal polygon of with 'width', that is at index 'i', 'j' in
-- a hexagonal tiling.
-- Called with arguments i, j, width
CREATE OR REPLACE FUNCTION hexagon_tile(INTEGER, INTEGER, FLOAT)
RETURNS GEOMETRY AS $$
    SELECT ST_TransScale(ST_GeomFromEWKT('POLYGON ((
         0.0   0.577350269189626, 
         0.5   0.288675134594813, 
         0.5  -0.288675134594813, 
         0.0  -0.577350269189626,
        -0.5  -0.288675134594813, 
        -0.5   0.288675134594813, 
         0.0   0.577350269189626 
        ))'),
        $1 + $2/2.0, $2*0.866025403784439, $3, $3)
$$ LANGUAGE SQL IMMUTABLE;


-- Returns the hexagon tile of width 'width' containing the cartesian point 'x', 'y'
-- Called with arguments x, y, width
CREATE OR REPLACE FUNCTION get_hexagon(FLOAT, FLOAT, FLOAT)
RETURNS GEOMETRY AS $$
    SELECT hexagon_tile(hexagon_i($1, $2, $3), hexagon_j($1, $2, $3), $3)
$$ LANGUAGE SQL IMMUTABLE;

-- Test the functions for generating the tiles
-- When exported using pgsql2shp, _test_hexagons should 
-- be a small patch of tiles, three across, and five high.
DROP TABLE IF EXISTS _test_hexagons;
CREATE TABLE _test_hexagons (
    i INTEGER,
    j INTEGER,
    hex GEOMETRY
);
INSERT INTO _test_hexagons VALUES  ( 0,-2, hexagon_tile( 0, -2, 15));
INSERT INTO _test_hexagons VALUES  ( 0,-1, hexagon_tile( 0, -1, 15));
INSERT INTO _test_hexagons VALUES  ( 0, 0, hexagon_tile( 0,  0, 15));
INSERT INTO _test_hexagons VALUES  ( 0, 1, hexagon_tile( 0,  1, 15));
INSERT INTO _test_hexagons VALUES  ( 0, 2, hexagon_tile( 0,  2, 15));
INSERT INTO _test_hexagons VALUES  (-1,-2, hexagon_tile(-1, -2, 15));
INSERT INTO _test_hexagons VALUES  (-1,-1, hexagon_tile(-1, -1, 15));
INSERT INTO _test_hexagons VALUES  (-1,-1, hexagon_tile(-1,  0, 15));
INSERT INTO _test_hexagons VALUES  (-1, 1, hexagon_tile(-1,  1, 15));
INSERT INTO _test_hexagons VALUES  (-1, 2, hexagon_tile(-1,  2, 15));
INSERT INTO _test_hexagons VALUES  ( 1,-2, hexagon_tile( 1, -2, 15));
INSERT INTO _test_hexagons VALUES  ( 1,-1, hexagon_tile( 1, -1, 15));
INSERT INTO _test_hexagons VALUES  ( 1,-1, hexagon_tile( 1,  0, 15));
INSERT INTO _test_hexagons VALUES  ( 1, 1, hexagon_tile( 1,  1, 15));
INSERT INTO _test_hexagons VALUES  ( 1, 2, hexagon_tile( 1,  2, 15));