Last active
October 27, 2024 07:01
-
-
Save mjumbewu/1761802ea06fb78c596f9cf8c9b2e769 to your computer and use it in GitHub Desktop.
PostGIS function to generate a grid of hexagonal cells in a PostgreSQL database
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
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; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment