sql - ST_HexagonGrid geom vector to find all points

Question

I was reviewing this function in PostGis

https://postgis.net/docs/manual-dev/ST_HexagonGrid.html

1) What I don't understand is what the underlying geom data would be. What's the source to get the USA map as shown? What is the DB schema? I think it could be one record, if I only need the USA boundary, not each state?

2) Is the result a list of points? or geom vectors?

3) If geom vectors, how do you convert them into points of lat and lng?

4) How to do approximate the hexigons to approximate a 50 mile radius from a point?

UPDATE:

I played with the width to try to get the correct number of hexagons based on Jim Jones example below. Unfortunately, something went wrong..

1) the length seems to have no relation to meters

2) There are multiple sized hexagons, which seems weird.

postgis_test=# WITH j AS (
postgis_test(# SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
postgis_test(#   SELECT 
postgis_test(#   generate_hexgrid(
postgis_test(#     5909968.8,
postgis_test(#     ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
postgis_test(# FROM usa_states)i) 
postgis_test-# SELECT count(j.hex) FROM j,usa_states
postgis_test-# WHERE ST_Intersects(usa_states.geom,j.hex);
 count 
-------
   119
(1 row)

postgis_test=# WITH j AS (
postgis_test(# SELECT ST_Transform((hex).geom,4326) AS hex FROM ( 
postgis_test(#   SELECT 
postgis_test(#   generate_hexgrid(
postgis_test(#     5909968.8,
postgis_test(#     ST_XMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMin(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_XMax(ST_Extent(ST_Transform(geom,3857))) ,
postgis_test(#     ST_YMax(ST_Extent(ST_Transform(geom,3857))) ) AS hex
postgis_test(# FROM usa_states)i) 
postgis_test-# SELECT DISTINCT st_area(j.hex) FROM j,usa_states
postgis_test-# WHERE ST_Intersects(usa_states.geom,j.hex);
     st_area      
------------------
 1219.78281686003
 2089.11341619338
 2089.11341619338
 3379.93344444246
  7051.4650344734
 12076.9943663072
(6 rows)

Answer 1

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:

• ST_Extent
• ST_Intersects