Hex grid algorithm for PostGIS

Hex binning is a great way to visualize data which is provided as polygon geometries. This article provides a simple PL/pgSQL script to generate a hexagonal grid from input geometries. A second query calculates values based on the intersections of grid cells and input geometries.

This is the first part of a workshop series which I hosted at FOSS4G 2017 in Boston. To all participants of that workshop. Sorry for the late post. I was meant to do this last month but got stuck in a post conference work stack.

The initial idea for this script comes from a method which is explained in detail in this blog post on The Academic Commons of The City University New York. I also received inspiration from this script by Rex Douglas. I am also aware of Cartos’ CDB_HexagonGrid() function, Postgis scripts by minus34 and the Hexagonal Grid function on the late DimensionalEdge blog. The last three functions I have not tested myself.

Requirements to my script were:

Cells generated on grids at different scales align in a regular pattern.

Generated cells must intersect the input geometries.

There should be no duplicate cells for neighbouring input geometries.

For this article I use Walhkreise of my home state of North Rhine-Westphalia with electoral results from the 2017 Bundestagswahl (General Election).

Wahlkreise of North Rhine-Westphalia

Here is the script with input parameters for a 10k regular hex grid.

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.

The generated hex grid overlaid to the input geometries looks like this.

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.

This process is designed for convenience not performance. It is easily doable to generate a stored function from this algorithm. The process slows down with huge series and it is recommended to use small input geometries or larger hex cells. I also run the process with static instead of temporary tables in order to improve performance. I have regularly generated grids with hundreds of millions of grid cells like this. Here a snapshot of an application showing population grid cells at a height of 500 meter per cell. You do the maths.

Population density grid in central Tokyo.

The next query generates a new table from the hex grid in which I calculate values by summing the area of the cell intersection divided by the area of the intersecting polygon times a selected field value from the input geometry. I group the sums by the indexed geometry field of the hex cells.

The performance of this query can be improved by updating fields in the existing table and calculating the intersection ratio separately. Here is what the distribution of voters looks like in NRW.

The regular pattern for a 5km grid overlaid to the 10km grid looks like this. This pattern is a requirement for dynamic web map layers.

As with all open source bits. We are more than happy if you use our algorithm. Please give credit to @GEOLYTIX as the author if you use our code in your work. Keep it free/libre.