Create hex grids over a layer

Question

I have a shapefile that has the outlines of the United States and I want to create hex grids that over that layer.

I tried using this function that will create hexes over a layer. I created an empty table and then added a geom column and then ran the custom function as shown.

CREATE TABLE IF NOT EXISTS hex_grid (hid serial NOT NULL PRIMARY KEY)
SELECT AddGeometryColumn('hex_grid', 'geom', 0, 'POLYGON', 2);

-- creating the hex grid
    SELECT genhexagons(0.05,ST_xmin(ST_Extent(SELECT geom FROM usa)),
    ST_ymin(ST_Extent(SELECT geom FROM usa)),
     ST_xmax(ST_Extent(SELECT geom FROM usa)),
     ST_ymax(ST_Extent(SELECT geom FROM usa)));

However I ran into this error that I'm not sure how to deal, I'm sure it deals with the creation of the hexes when it exceeds the US outline boundary. See below:

Answer 1

Looking at the post you linked to, the genhexagons function takes 5 float parameters, cell width, and the (x,y) coords for SW and NE corners of extent. You're passing in a float and a geometry.

CREATE OR REPLACE FUNCTION genhexagons(width float, xmin float,ymin  float,xmax float,ymax float)

You might be able to use ST_xmin(ST_Extent(your_geom)) (and ST_ymin, ST_ymax etc.) if you don't want to hard-code the coordinate values of the extent in your function call.

I've tried the code in that post verbatim, and it works as advertised - it creates a hex grid over Kenya.

EDIT

Might be simpler if you specify the bounds manually. For the continguous US, I'm guessing the corners of the extent are approximately (-124.8,25.10) and (-66.9,49.6)

so calling

genhexagons(1.0, -124.8, 25.10, -66.9, 49.6)

gives 1 degree grid over the contiguous US...

However it appears this code will generate too many rows of tiles... I find it covers most of Canada too :/

I found I was able to clip the results to the bounding box by dropping any tiles which didn't overlap the bounding box, and writing these to a new table called clipped.

create table clipped as (
    select * from hex_grid where 
        ST_Intersects(
            st_setSRID(the_geom,4326),
            st_setSRID(st_envelope('POLYGON((-124.8 25.10, -124.8 49.6, -66.9 49.6, -66.9 25.10, -124.8 25.10))'::geometry),4326)
    )
);

Bringing this into QGIS...

Answer 2

I wrote a PostGIS function to generate hex grids on top of another layer.

DO $$
DECLARE
  _curs   CURSOR FOR SELECT geom3857 FROM nrw;
  _table  TEXT    := 'nrw_hx_10k';
  _srid   INTEGER := 3857;
  _height NUMERIC := 10000;
  _width  NUMERIC := _height * 0.866;
  _geom   GEOMETRY;
  _hx     TEXT    := 'POLYGON((' || 0 || ' ' || 0 || ',' || (_width * 0.5) || ' ' || (_height * 0.25) || ',' ||
                 (_width * 0.5) || ' '
                 || (_height * 0.75) || ',' || 0 || ' ' || _height || ',' || (-1 * (_width * 0.5)) || ' ' ||
                 (_height * 0.75) || ',' ||
                 (-1 * (_width * 0.5)) || ' ' || (_height * 0.25) || ',' || 0 || ' ' || 0 || '))';
  _hx_g   GEOMETRY := ST_SetSRID(_hx::GEOMETRY, _srid);

BEGIN
  CREATE TEMP TABLE hx_tmp (geom GEOMETRY(POLYGON));

  OPEN _curs;
  LOOP
    FETCH
    _curs INTO _geom;
    EXIT WHEN NOT FOUND;

    INSERT INTO hx_tmp
      SELECT
        ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON) geom
      FROM
        generate_series(
          (st_xmin(_geom) / _width)::INTEGER * _width - _width,
          (st_xmax(_geom) / _width)::INTEGER * _width + _width,
          _width) x_series,
        generate_series(
          (st_ymin(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) - _height,
          (st_ymax(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) + _height,
          _height * 1.5) y_series
      WHERE
        ST_Intersects(ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON), _geom);

    INSERT INTO hx_tmp
      SELECT ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON) geom
      FROM
        generate_series(
          (st_xmin(_geom) / _width)::INTEGER * _width - (_width * 1.5),
          (st_xmax(_geom) / _width)::INTEGER * _width + _width,
          _width) x_series,
        generate_series(
          (st_ymin(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) - (_height * 1.75),
          (st_ymax(_geom) / (_height * 1.5))::INTEGER * (_height * 1.5) + _height,
          _height * 1.5) y_series
      WHERE
        ST_Intersects(ST_Translate(_hx_g, x_series, y_series)::GEOMETRY(POLYGON), _geom);

  END LOOP;
  CLOSE _curs;

  CREATE INDEX sidx_hx_tmp_geom ON hx_tmp USING GIST (geom);
  EXECUTE 'DROP TABLE IF EXISTS '|| _table;
  EXECUTE 'CREATE TABLE '|| _table ||' (geom GEOMETRY(POLYGON, '|| _srid ||'))';
  EXECUTE 'INSERT INTO '|| _table ||' SELECT * FROM hx_tmp GROUP BY geom';
  EXECUTE 'CREATE INDEX sidx_'|| _table ||'_geom ON '|| _table ||' USING GIST (geom)';
  DROP TABLE IF EXISTS hx_tmp;
END $$;

Input parameters must be set for: _curs: The geometry field name and the table name of input geometries. _table: The name of the output table. _srid: The geographic projection code of the input (and output) geometries. _height: The height of a hexgon grid cell in projection units.

I generate a scaled hexagon geometry in the declare block and then loop through the input geometries. In the loop, I generate series for the x and y extent plus some for each input geometry. The hexagon is translated and inserted into a temporary table if the two geometries intersect. A second pair of series is generated alternative offset rows. Finally I group the hexagon grid cells by their geometries to remove duplicates.

There is a more detailed description and some background on medium.

Answer 3

A simple Function for a hex grid Github Repo: https://github.com/imran-5/Postgis-Custom

Function====================================================================

CREATE OR REPLACE FUNCTION public.I_Grid_Hex(geom geometry, radius double precision)
    RETURNS SETOF geometry AS $BODY$
    DECLARE
    srid INTEGER := 4326;
    input_srid INTEGER;
    x_max DECIMAL;
    y_max DECIMAL;
    x_min DECIMAL;
    y_min DECIMAL;
    x_series DECIMAL;
    y_series DECIMAL;
    b float :=radius/2;
    a float :=b/2; --sin(30)=.5
    c float :=2*a;
--temp     GEOMETRY := ST_GeomFromText(FORMAT('POLYGON((0 0, %s %s, %s %s, %s %s, %s %s, %s %s, 0 0))',
--                       (b), (a), (b), (a+c), (0), (a+c+a), (-1*b), (a+c), (-1*b), (a)), srid);
geom_grid     GEOMETRY := ST_GeomFromText(FORMAT('POLYGON((0 0, %s %s, %s %s, %s %s, %s %s, %s %s, 0 0))',
                      (radius *  0.5), (radius * 0.25),
                      (radius *  0.5), (radius * 0.75),
                                   0 ,  radius,
                      (radius * -0.5), (radius * 0.75),
                      (radius * -0.5), (radius * 0.25)), srid);
BEGIN
CASE st_srid(geom) WHEN 0 THEN
    geom := ST_SetSRID(geom, 4326);
    RAISE NOTICE 'SRID Not Found.';
ELSE
    RAISE NOTICE 'SRID Found.';
END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);
    x_series := ceil ( @( x_max - x_min ) / radius);
    y_series := ceil ( @( y_max - y_min ) / radius);
RETURN QUERY
    With foo as(SELECT
        ST_Translate ( cell, x*(2*a+c)+x_min, y*(2*(c+a))+y_min) AS hexa
    FROM
        generate_series ( 0, x_series, 1) AS x,
        generate_series ( -1, y_series, 1) AS y,
        (
        SELECT ST_Translate(geom_grid::geometry, b , a+c)  as cell
            union
        SELECT geom_grid AS cell
        ) AS foo
    ) select ST_CollectionExtract(ST_Collect(ST_Transform(ST_Intersection(ST_CollectionExtract(hexa, 3), geom), input_srid)),3)
    from foo where ST_Intersects(hexa, geom);
END;
$BODY$ LANGUAGE 'plpgsql' VOLATILE;

Query for above question:

SELECT I_Grid_Hex(st_envelope(st_collect(geom)), 1) from "us-states"

Result

Another simple query

SELECT I_Grid_Hex(geom, .0001) from polygons limit 1

Result======================================================================