create-hex-grid.sql

-- This script is available for use under the BSD License
-- 
-- Jonas Häggqvist <[email protected]>
-- Guillaume Rischard <[email protected]>
-- Create grid, based on the function from http://trac.osgeo.org/postgis/wiki/UsersWikiGenerateHexagonalGrid
-- Hexagons are cut using the country polygon as a cookie cutter.

-- Magic values:

-- 2169 is the SRID.
-- 48000, 107000, 57000, 139000 are the bbox for the country in EPSG 2169.
-- 35040 is the negative x offset for the bbox when we multiply the x value by sqrt(3).

-- drop previous grid
DROP TABLE IF EXISTS "grid";

SELECT * INTO grid
    FROM
    (
        SELECT
            ST_Intersection(fullgrid.geometry, country.way) AS geometry
        FROM (
            select way from planet_osm_polygon WHERE admin_level='2' ) as country,
            (
                WITH    
                cell as (
                    select
                        1000/2 as size, -- Radius of the hexagon (which is then divided by 2)
                        1.73205080757 as xyratio -- horizontal spacing = vertical spacing * sqrt(3)
                ),
                hexagon as (
                    SELECT ST_GeomFromText(
                            'POLYGON((
                                0 0,
                                '||(select size*xyratio from cell)||' '||(select size from cell)||',
                                '||(select size*xyratio from cell)||' '||(select size*3 from cell)||',
                                '||'0 '||(select size*4 from cell)||',
                                -'||(select size*xyratio from cell)||' '||(select size*3 from cell)||',
                                -'||(select size*xyratio from cell)||' '||(select size from cell)||',
                                0 0))'
                            ) as the_geom
                        )
                SELECT
                    ST_SetSRID(st_translate(the_geom, x_series*(select xyratio from cell)-35040, y_series), 2169) as geometry
                FROM
                    generate_series(48000, 107000, (select size*2 from cell)) as x_series,
                    generate_series(57000, 139000, (select size*6 from cell)) as y_series,
                    (
                        SELECT the_geom from hexagon
                        UNION
                        SELECT ST_Translate((SELECT the_geom from hexagon), (select size*xyratio from cell), (select size*3 from cell))as the_geom
                    ) as two_hex
            ) AS fullgrid
        ) as geometry;