How to create hexagonal grids in PostGIS

gistfile1.txt

-- Function: makegrid_epsg3857(text, text, text, numeric)

-- DROP FUNCTION makegrid_epsg3857(text, text, text, numeric);

CREATE OR REPLACE FUNCTION makegrid_epsg3857(schemaname text, boundingbox text, gridtable text, halfwidth numeric)
  RETURNS text AS
$BODY$
    DECLARE
    tbl_cnt int;
    XMIN numeric;
    XMAX numeric;
    YMIN numeric;
    YMAX numeric;
    x_value numeric;
    y_value numeric;
    x_count integer;
    y_count integer;
    y_offset numeric;
    y_value_adj numeric;
    sequencevar text := gridtable ||'_ogc_fid_seq';

    BEGIN
    -- Check to see if grid table already exists
    SELECT COUNT(*) INTO tbl_cnt FROM information_schema.tables WHERE table_schema = schemaname AND table_name = gridtable;
    -- If grid table already exists, drop table
    IF (tbl_cnt > 0) THEN
    EXECUTE 'DROP TABLE ' || gridtable;
    ELSE
    END IF;

    -- Create grid table
    EXECUTE 'CREATE TABLE ' || gridtable ||'
    (
    ogc_fid serial,
    x_count integer,
    x numeric(8,1),
    y_count integer,
    y numeric(8,1),
    geometry geometry(Polygon,3857)
    )
    WITH (
      OIDS=FALSE
    );

    CREATE INDEX ' || gridtable || '_geom
     ON ' || gridtable ||'
     USING gist
     (geometry);';

    -- load extents of the bounding box and values for the initial y parameters
    EXECUTE 'SELECT floor(st_xmin(ST_transform(geometry,3857))) FROM '|| boundingbox INTO XMIN;
    EXECUTE 'SELECT ceiling(st_xmax(ST_transform(geometry,3857))) FROM '|| boundingbox INTO XMAX;
    EXECUTE 'SELECT floor(st_ymin(ST_transform(geometry,3857))) FROM '|| boundingbox INTO YMIN;
    EXECUTE 'SELECT ceiling(st_ymax(ST_transform(geometry,3857))) FROM '|| boundingbox INTO YMAX;
    y_count = 1;
    y_value = YMIN;

    LOOP
    -- for each y value, reset x to XMIN and subloop through the x values
        x_count = 1;
        x_value = XMIN;
        LOOP
        -- for every even numbered x_count, offset y by half the height of the hexagon
        y_offset = (ceiling((x_count+1)/2::numeric) - floor((x_count+1)/2::numeric))*round(halfwidth*SQRT(3)::numeric,1)/2;
        -- add the offset to the y_value
        y_value_adj = y_value + y_offset;
        EXECUTE 'INSERT INTO ' || gridtable || ' VALUES(NEXTVAL('|| quote_literal(sequencevar) ||'), '|| x_count ||', '|| x_value ||', ' || y_count ||', ' || y_value_adj ||', NULL)';
        x_count = x_count + 1;
        x_value = x_value + round(halfwidth*3/2::numeric,1);
        EXIT WHEN x_value > XMAX;
        END LOOP;
    -- after exiting the subloop, increment the y count and y value
    y_count = y_count + 1;
    y_value = y_value + round(halfwidth*SQRT(3)::numeric,1);
    EXIT WHEN y_value > YMAX;
    END LOOP;

    -- With the table now populated with x y points, the last step is to create a hexagon geometry for each.
    EXECUTE 'UPDATE '|| gridtable ||' SET geometry = ST_transform(ST_Polygon(ST_makeline(ARRAY[st_makepoint(x-(0.5*'|| halfwidth ||'), y+(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x-(1*'|| halfwidth ||'), y), st_makepoint(x-(0.5*'|| halfwidth ||'), y-(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x+(0.5*'|| halfwidth ||'), y-(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x+(1*'|| halfwidth ||'), y), st_makepoint(x+(0.5*'|| halfwidth ||'), y+(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x-(0.5*'|| halfwidth ||'), y+(SQRT(3)*0.5*'|| halfwidth ||'))]),3857),3857)';

    RETURN 'HEXAGON GRID CREATED';
    END;
    $BODY$
  LANGUAGE plpgsql VOLATILE STRICT
  COST 100;

-- USAGE: 
-- SELECT makegrid_epsg3857('public',E'(SELECT ST_SetSRID(ST_MakePolygon(ST_GeomFromText(\'LINESTRING(-19.6875 26.1949,5.7129 26.1949,5.7129 44.2137,-19.6875 44.2137,-19.6875 26.1949)\')),4326) as geometry) a', 'hexamap', 10750)

-- REFERENCES:
-- http://www.dimensionaledge.com/main/postgis/how-to-create-hexagonal-grids-in-postgis/