According to the author, the following function should create a grid with the extent based on the given BBOX and the cell size in meters.
SRID 3857 units are [very approximately] meters, and using this
projection will create hex cells that "look right" on a web map (most
of which use a web mercator projection).
CREATE OR REPLACE FUNCTION generate_hexgrid(width float, xmin float, ymin float, xmax float, ymax float, srid int default 3857)
RETURNS TABLE(gid text,geom geometry(Polygon)) AS $$
DECLARE
b float := width / 2;
a float := tan(radians(30)) * b;
c float := 2 * a;
height float := 2 * (a + c);
index_xmin int := floor(xmin / width);
index_ymin int := floor(ymin / height);
index_xmax int := ceil(xmax / width);
index_ymax int := ceil(ymax / height);
snap_xmin float := index_xmin * width;
snap_ymin float := index_ymin * height;
snap_xmax float := index_xmax * width;
snap_ymax float := index_ymax * height;
ncol int := abs(index_xmax - index_xmin);
nrow int := abs(index_ymax - index_ymin);
polygon_string varchar :=
'POLYGON((' || 0 || ' ' || 0 || ' , ' || b || ' ' || a || ' , ' ||
b || ' ' || a + c || ' , ' || 0 || ' ' || a + c + a || ' , ' ||
-1 * b || ' ' || a + c || ' , ' || -1 * b || ' ' || a || ' , ' ||
0 || ' ' || 0 ||'))';
BEGIN
RETURN QUERY
SELECT
format('%s %s %s', width,
x_offset + (1 * x_series + index_xmin),
y_offset + (2 * y_series + index_ymin)),
ST_SetSRID(ST_Translate(two_hex.geom,
x_series * width + snap_xmin,
y_series * height + snap_ymin), srid)
FROM generate_series(0, ncol, 1) AS x_series,
generate_series(0, nrow, 1) AS y_series,
(SELECT 0 AS x_offset, 0 AS y_offset, polygon_string::geometry AS geom
UNION
SELECT 0 AS x_offset, 1 AS y_offset, ST_Translate(polygon_string::geometry, b , a + c) AS geom
) AS two_hex;
END; $$ LANGUAGE plpgsql;
Considering you have a table called usa and it contains the geometries of this shapefile you should be able to create a grid that overlaps the USA map with the following query:
CREATE TABLE usa_hex AS
WITH j AS (
SELECT ST_Transform((hex).geom,4326) AS hex FROM (
SELECT
generate_hexgrid(
80467,
ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
FROM usa)i)
SELECT j.hex FROM j,usa
WHERE ST_Intersects(usa.geom,j.hex);
EDIT: It is still not an answer, as it does not create the hexagons using meters, but it might help other users. The following function (derived from this answer) creates geometry type hexagons of exact the same size in degrees.
CREATE OR REPLACE FUNCTION generate_hexagons(width FLOAT, bbox BOX2D, srid INTEGER DEFAULT 4326)
RETURNS TABLE (gid INTEGER, hexagon GEOMETRY) AS $$
DECLARE
b FLOAT := width/2;
a FLOAT := b/2;
c FLOAT := 2*a;
height FLOAT := 2*a+c;
ncol FLOAT := ceil(abs(ST_Xmax(bbox)-ST_Xmin(bbox))/width);
nrow FLOAT := ceil(abs(ST_Ymax(bbox)-ST_Ymin(bbox))/height);
polygon_string VARCHAR := 'POLYGON((' ||
0 || ' ' || 0 || ' , ' || b || ' ' || a || ' , ' || b || ' ' || a+c || ' , ' || 0 || ' ' || a+c+a || ' , ' ||
-1*b || ' ' || a+c || ' , ' || -1*b || ' ' || a || ' , ' || 0 || ' ' || 0 || '))';
BEGIN
RETURN QUERY
SELECT
row_number() OVER ()::INTEGER,
ST_SetSRID(
ST_Translate(geom, x_series*(2*a+c)+ST_Xmin(bbox), y_series*(2*(c+a))+ST_Ymin(bbox)),srid)
FROM generate_series(0, ncol::INTEGER, 1) AS x_series,
generate_series(0, nrow::INTEGER,1 ) AS y_series,
(SELECT polygon_string::GEOMETRY AS geom
UNION
SELECT ST_Translate(polygon_string::GEOMETRY, b, a + c) AS geom) AS two_hex;
END;
$$ LANGUAGE plpgsql;
Overlapping with the data set used above:
WITH j (hex_rec) AS (
SELECT generate_hexagons(3.0,ST_Extent(geom))
FROM usa
)
SELECT (hex_rec).gid,(hex_rec).hexagon FROM j, usa
WHERE ST_Intersects(usa.geom,(hex_rec).hexagon);
Further reading:
