Spatial Hexagon Binning in POSTGIS

FacebookTwitterGoogle+RedditLinkedIn

hex4Hexagon binning is an important tool for visualizing spatial relationships but it’s a pain to implement. Current solutions involve partial implementations for PostGIS and work arounds based on generating hexagons in QGIS with the MMQGIS plugin. In this brief post, I provide a SQL function for quickly generating hexagon polygon layers in PostGIS of any size, region of interest, and projection.

 

The PL/pgSQL  function below takes in 5 parameters, the width of a regular hexagon and four parameters providing the desired geographic extent. It is a heavily modified version of the code provided here, and works by constructing two contiguous hexagons and then tiling them enough times to cover the entire region of interest. By default it places the results in a new table “hex_grid” and uses the World Geodetic System 84 (SRID 4326). The unit of measurement is then in degrees latitude and longitude.

Code

CREATE TABLE hex_grid (gid serial not null primary key);
SELECT addgeometrycolumn('hex_grid','the_geom', 0, 'POLYGON', 2); 

CREATE OR REPLACE FUNCTION genhexagons(width float, xmin float,ymin  float,xmax float,ymax float  )
RETURNS float AS $total$
declare
	b float :=width/2;
	a float :=b/2; --sin(30)=.5
	c float :=2*a;
	height float := 2*a+c;  --1.1547*width;
	ncol float :=ceil(abs(xmax-xmin)/width);
	nrow float :=ceil(abs(ymax-ymin)/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
    INSERT INTO hex_grid (the_geom) SELECT st_translate(the_geom, x_series*(2*a+c)+xmin, y_series*(2*(c+a))+ymin)
    from generate_series(0, ncol::int , 1) as x_series,
    generate_series(0, nrow::int,1 ) as y_series,
    (
       SELECT polygon_string::geometry as the_geom
       UNION
       SELECT ST_Translate(polygon_string::geometry, b , a+c)  as the_geom
    ) as two_hex;
    ALTER TABLE hex_grid
	ALTER COLUMN the_geom TYPE geometry(Polygon, 4326)
	USING ST_SetSRID(the_geom,4326);
    RETURN NULL;
END;
$total$ LANGUAGE plpgsql;

--width in the units of the projection, xmin,ymin,xmax,ymax
SELECT genhexagons(0.1,34.0,-5.0,42.0,1.25);

Example

The final product is a PostGIS polygon layer of contiguous hexagons covering the region of interest. In an example below, I’ve created a polygon layer with hexagons 11 km in width. I imported that layer into QGIS along with road coverage for Kenya. Using the segment length tool, I tally the the total length of roads criss-crossing each hexagon. In this application, the hexagon binning nicely captures a coarse estimate of local infrastructure that can be easily combined with other spatial data.

This slideshow requires JavaScript.

FacebookTwitterGoogle+RedditLinkedIn

6 comments to Spatial Hexagon Binning in POSTGIS

  • […] For these examples, I’m using a base layer consisting of a hexagon lattice covering the 1950s borders of Kenya.  It consists of 6,300 regular hexagons each 11km wide.  I created this lattice and added it as a table called “hex_grid” using a function described in my previous post. […]

  • andrew lainton

    Hi Rex I modified your script (my sql is terrible) to use the Eckhert IV projection SRID=4326 as follows and it gives an error. What we needed to do was create hexagons over a very large area which needed to be exactly the same in area to avoid modelling distortions and so for that reason create them in a spherical equal area projection, in this case eaxctly 1/7 sqm in area over the bounding box of the British national grid.

    I cant see the error but my experience in postgis sql is minimal

    CREATE TABLE hex_grid (gid serial not null primary key);
    SELECT addgeometrycolumn('hex_grid','the_geom', 4326, 'POLYGON', 2);
    CREATE OR REPLACE FUNCTION genhexagons(width float, xmin float,ymin float,xmax float,ymax float )
    RETURNS float AS $total$
    declare
    b float :=width/2;
    a float :=b/2; --sin(30)=.5
    c float :=2*a;
    height float := 2*a+c; --1.154700538379252*width;
    ncol float :=ceil(abs(xmax-xmin)/width);
    nrow float :=ceil(abs(ymax-ymin)/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
    INSERT INTO hex_grid (the_geom) SELECT st_translate(the_geom, x_series*(2*a+c)+xmin, y_series*(2*(c+a))+ymin)
    from generate_series(0, ncol::int , 1) as x_series,
    generate_series(0, nrow::int,1 ) as y_series,
    (
    SELECT polygon_string::geometry as the_geom
    UNION
    SELECT ST_Translate(polygon_string::geometry, b , a+c) as the_geom
    ) as two_hex;
    ALTER TABLE hex_grid
    ALTER COLUMN the_geom TYPE geometry(Polygon, 4326)
    USING ST_SetSRID(the_geom,4326);
    RETURN NULL;
    END;
    $total$ LANGUAGE plpgsql;
    --width in the units of the projection, xmin,ymin,xmax,ymax
    SELECT genhexagons(406.15,-661595.4844,6041382.8902,195092.365400001,7119317.9239);

    ERROR: Geometry SRID (0) does not match column SRID (4326)
    SQL state: 22023
    Context: SQL statement "INSERT INTO hex_grid (the_geom) SELECT st_translate(the_geom, x_series*(2*a+c)+xmin, y_series*(2*(c+a))+ymin)
    from generate_series(0, ncol::int , 1) as x_series,
    generate_series(0, nrow::int,1 ) as y_series,
    (
    SELECT polygon_string::geometry as the_geom
    UNION
    SELECT ST_Translate(polygon_string::geometry, b , a+c) as the_geom
    ) as two_hex"
    PL/pgSQL function genhexagons(double precision,double precision,double precision,double precision,double precision) line 19 at SQL statement

    That would make me very happy

    Also - dont worry your time if your busy- would like to create a,b,c coordinates (s, p and c in our project) as described here under

    http://www.redblobgames.com/grids/hexagons/

    Cube coordinates

    As this would make the maths for calculating coordinates and distances in the project very fast.

    Finally, though could be done outside, would be super nice to generate points at the same time at the centroids with the unique ID and s,p,c coordinates, to join to, as labelling and managing these would be so much easier than with geometries of 5 million + polygons.

  • Hi Rex, I like your code. This site was a great help to write my own function. Something however doesn't stick right. The width of a hex should be 3/4 the height. Making a = height * 1/4, a + c = height * 3/4, and a+c+a = height. Similar b can be expressed as width/2 thus getting completely rid of these auxiliary terms.

    • David

      Edited the code to make it work for equal side lengths where you only have to enter the side length. Used the theory from http://www.redblobgames.com/grids/hexagons/ to complete.

      DROP TABLE IF EXISTS hex_grid_;
      CREATE TABLE hex_grid (
      gid serial not null primary key,
      the_geom geometry(polygon, 3400) NOT NULL
      );

      CREATE INDEX gix_hex_grid ON hex_grid USING GIST(the_geom);
      SELECT genhexagons(100,
      (select st_xmin(geom) from (select st_extent(geom) geom from table_X)foo),
      (select st_ymin(geom) from (select st_extent(geom) geom from table_X)foo),
      (select st_xmax(geom) from (select st_extent(geom) geom from table_X)foo),
      (select st_ymax(geom) from (select st_extent(geom) geom from table_X)foo));

      CREATE OR REPLACE FUNCTION genhexagons (size float, xmin float,ymin float,xmax float,ymax float)
      RETURNS float AS $total$
      declare

      sx1 float := cos((60 * 1 + 30) * (pi()/180)) * size;
      sx2 float := cos((60 * 2 + 30) * (pi()/180)) * size;
      sx3 float := cos((60 * 3 + 30) * (pi()/180)) * size;
      sx4 float := cos((60 * 4 + 30) * (pi()/180)) * size;
      sx5 float := cos((60 * 5 + 30) * (pi()/180)) * size;
      sx6 float := cos((60 * 6 + 30) * (pi()/180)) * size;

      sy1 float := sin((60 * 1 + 30) * (pi()/180)) * size;
      sy2 float := sin((60 * 2 + 30) * (pi()/180)) * size;
      sy3 float := sin((60 * 3 + 30) * (pi()/180)) * size;
      sy4 float := sin((60 * 4 + 30) * (pi()/180)) * size;
      sy5 float := sin((60 * 5 + 30) * (pi()/180)) * size;
      sy6 float := sin((60 * 6 + 30) * (pi()/180)) * size;

      width float := ((size / 2) * sqrt(3)) * 2;
      a float := cos((60) * (pi()/180)) * size;
      b float := width / 2;
      c float := 2 * a;

      height float := size + c;

      ncol float :=ceil(abs(xmax-xmin)/width);

      nrow float :=ceil(abs(ymax-ymin)/height);

      polygon_string varchar := 'SRID=3400;POLYGON((' ||

      sx1||' '||sy1||' , '||
      sx2||' '||sy2||' , '||
      sx3||' '||sy3||' , '||
      sx4||' '||sy4||' , '||
      sx5||' '||sy5||' , '||
      sx6||' '||sy6||' , '||
      sx1||' '||sy1||
      '))';

      BEGIN
      INSERT INTO hex_grid (the_geom) SELECT st_translate(the_geom, x_series*(width)+xmin, y_series*(2*(c+a))+ymin)
      from generate_series(0, ncol::int,1) as x_series,
      generate_series(0, nrow::int,1) as y_series,
      (
      SELECT polygon_string::geometry as the_geom

      UNION
      SELECT ST_Translate(polygon_string::geometry, b , a+c) as the_geom
      ) as two_hex;

      RETURN NULL;
      END;
      $total$ LANGUAGE plpgsql

  • Angus Carr

    I modified your script to generate hexagons of a given area. I noticed that the geometry of your hexagons was a bit off (They seem vertically squished) so I went back to the geometry of the hexagon to generate correct polygons. Fundamentally, I calculate the width from the area, then use your system. I wanted to check the area to make sure I got the width function correct, and found the area was 0.86 of what it should have been. Not a problem for your application, but I thought I should share.

    CREATE OR REPLACE FUNCTION genhexagons(area float, xmin float,ymin float,xmax float,ymax float )
    RETURNS table ( shape geometry(POLYGON) ) AS
    $total$
    DECLARE
    -- Geometry symbols and formulae based on http://www.calculatorsoup.com/calculators/geometry-plane/polygon.php
    width float = ( ( area / tan (pi() / 6.0 ) / 6.0 ) ^ 0.5 ) * 2.0 ;
    b float =width/2.0; -- apothem (r) -- inner radius
    a float = tan(pi() / 6) * b ; -- R/2 = tan (pi/6) * r
    c float =2.0*a; -- R --Outer radius, also perimeter /6
    height float := 2.0*a+c; --1.1547*width;
    ncol float :=ceil(abs(xmax-xmin)/width);
    nrow float :=ceil(abs(ymax-ymin)/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 st_translate(the_geom, x_series*(width)+xmin, y_series*(2*(c+a))+ymin)
    from generate_series(0, ncol::int , 1) as x_series,
    generate_series(0, nrow::int,1 ) as y_series,
    (
    SELECT polygon_string::geometry as the_geom
    UNION
    SELECT ST_Translate(polygon_string::geometry, b , a+c) as the_geom
    ) as two_hex;

    END;
    $total$
    LANGUAGE plpgsql;
    --width in the units of the projection, xmin,ymin,xmax,ymax

    --DELETE FROM hex_grid ;
    --CREATE TABLE hex_grid (gid serial not null primary key, area float, shape geometry(POLYGON,3161));
    INSERT INTO hex_grid (area, shape )
    SELECT ST_Area(shape),ST_SetSRID(shape,3161) FROM genhexagons(10000.0,34.0,-5.0,42.0,1.25) ;

  • I did a completely rewrite... Generating a series with the correct centroids was much faster than applying the shift to the series value. I also use a test table to get the extent as well as to check whether a grid cell is intersecting the test geometry. I had a problem where a complete grid was way to big for my server to handle. It's important that the test geometry is indexed. I also prefer to run a SQL query instead of creating a function first and then run the function. I got rid of the union and created the alt lines separate in order to check the individual hex cells for intersection. I multiply with decimals instead dividing by fractions. Just something which makes things easier to read for myself. Here is the code. Might be helpful to some.

    SELECT ST_Translate(ST_SetSRID(geom,3857), x_series, y_series)::geometry(Polygon,3857) geom_3857
    INTO _table_out
    FROM
    generate_series((SELECT st_xmin(st_extent(st_transform(geom_4326, 3857)))::int FROM _table_in), (SELECT st_xmax(st_extent(st_transform(geom_4326, 3857)))::int FROM _table_in), _width) x_series,
    generate_series((SELECT st_ymin(st_extent(st_transform(geom_4326, 3857)))::int FROM _table_in), (SELECT st_ymax(st_extent(st_transform(geom_4326, 3857)))::int FROM _table_in), _height * 1.5) y_series,
    (SELECT ('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 || '))')::geometry geom) SQ,
    test_poly
    WHERE ST_Intersects(ST_Translate(ST_SetSRID(geom,3857), x_series, y_series)::geometry(Polygon,3857), _table_in._geom_in);

    INSERT INTO _table_out
    SELECT ST_Translate(ST_SetSRID(geom,3857), x_series, y_series)::geometry(Polygon,3857) geom_3857
    FROM
    generate_series((SELECT st_xmin(st_extent(st_transform(geom_4326, 3857)))::int FROM test_poly) - _width * 0.5, (SELECT st_xmax(st_extent(st_transform(geom_4326, 3857)))::int FROM test_poly), _width) x_series,
    generate_series((SELECT st_ymin(st_extent(st_transform(geom_4326, 3857)))::int FROM test_poly) - _height * 0.75, (SELECT st_ymax(st_extent(st_transform(geom_4326, 3857)))::int FROM test_poly), _height * 1.5) y_series,
    (SELECT ('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 || '))')::geometry geom) SQ,
    test_poly
    WHERE ST_Intersects(ST_Translate(ST_SetSRID(geom,3857), x_series, y_series)::geometry(Polygon,3857), _table_in._geom_in);

Leave a Reply to Dennis Cancel reply

  

  

  

You can use these HTML tags

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>