PostGIS function to generate a grid of hexagonal cells in a PostgreSQL database

1-hexgrid_illustration.svg

2-generate_hexgrid.sql

/*
The default SRID is EPSG 3857 (web mercator -- https://epsg.io/3857). However
you can use any SRID you want. All input parameters should be interpreted as
coordinates and distances in whatever the SRID is set to.

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).

If you have bounds in lat/lng degrees, you can convert those into web mercator.
To use EPSG 4326 (geodetic latitude and longitude -- https://epsg.io/4326)
degrees as the bounds, you can do the following:

    SELECT gid, ST_Transform(geom, 4326) AS geom
    FROM generate_hexgrid(

      -- Width of cell, in meters
      8192,

      -- Minimum x and y
      ST_X(ST_Transform(ST_SetSRID(ST_GeomFromText('POINT(-75.60447692871092 39.782685009007075)'), 4326), 3857)),
      ST_Y(ST_Transform(ST_SetSRID(ST_GeomFromText('POINT(-75.60447692871092 39.782685009007075)'), 4326), 3857)),

      -- Maximum x and y
      ST_X(ST_Transform(ST_SetSRID(ST_GeomFromText('POINT(-74.78736877441406 40.159459579477925)'), 4326), 3857)),
      ST_Y(ST_Transform(ST_SetSRID(ST_GeomFromText('POINT(-74.78736877441406 40.159459579477925)'), 4326), 3857)),

      -- The input SRID, default 3857
      3857
    );

The geometry returned from this function also uses EPSG 3857 coordinates, or
whatever the input SRID is, hence the use of an additional ST_Transform in the
SELECT above.

The gid should be unique for (and characteristic to) each cell. In other words,
If you run this function twice with two distinct but overlapping bounding boxes
using the same cell width, the cells that overlap should have the same gid. So,
if you INSERT these cells into a table with a unique gid column, you should be
able to ignore conflicts (ON CONFLICT DO NOTHING).

Adapted from http://rexdouglass.com/spatial-hexagon-binning-in-postgis/
Snapping inspired by https://medium.com/@goldrydigital/hex-grid-algorithm-for-postgis-4ac45f61d093
*/

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 $grid$
declare
  b float := width / 2;
  a float := tan(radians(30)) * b;  -- tan(30) = 0.577350269
  c float := 2 * a;

  -- NOTE: The height of one cell is (2a + c), or about 1.154700538 * width.
  --       however, for each row, we shift vertically by (2[a + c]) to properly
  --       tesselate the hexagons. Thus, to determine the number of rows needed,
  --       we use the latter formula as the height of a row.
  height float := 2 * (a + c);

  -- Snap the min/max coords to a global grid according to the cell width, so
  -- that we minimize the chances of generating misaligned grids for overlapping
  -- regions.
  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;

  -- Calculate the total number of columns and rows. Note that the number of
  -- rows is actually half the number of rows, since each vertical iteration
  -- accounts for two "rows".
  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

    -- gid is made of the global x index of the cell, the global y index of the
    -- cell, and the cell width.
    format('%s %s %s',
           width,
           x_offset + (1 * x_series + index_xmin),
           y_offset + (2 * y_series + index_ymin)),

    -- geom is transformed using the width and height of a series, and set to
    -- the SRID specified.
    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,

    -- two_hex is a pair of hex cells, one roughly below the other. Thus, both
    -- have an x_offset of 0, but the second has a y_offset of 1.
    (
      -- Series cell #1
      SELECT
        0 AS x_offset,
        0 AS y_offset,
        polygon_string::geometry AS geom

      UNION
      
       -- Series cell #2
      SELECT
        0 AS x_offset,
        1 AS y_offset,
        ST_Translate(polygon_string::geometry, b , a + c)  AS geom
    ) AS two_hex;
END;
$grid$ LANGUAGE plpgsql;