Skip to content

Instantly share code, notes, and snippets.

@vehrka
Created April 14, 2015 07:23
Show Gist options
  • Save vehrka/c33c65ab366e7aec14ad to your computer and use it in GitHub Desktop.
Save vehrka/c33c65ab366e7aec14ad to your computer and use it in GitHub Desktop.
How to create hexagonal grids in PostGIS
-- Function: makegrid_epsg3857(text, text, text, numeric)
-- DROP FUNCTION makegrid_epsg3857(text, text, text, numeric);
CREATE OR REPLACE FUNCTION makegrid_epsg3857(schemaname text, boundingbox text, gridtable text, halfwidth numeric)
RETURNS text AS
$BODY$
DECLARE
tbl_cnt int;
XMIN numeric;
XMAX numeric;
YMIN numeric;
YMAX numeric;
x_value numeric;
y_value numeric;
x_count integer;
y_count integer;
y_offset numeric;
y_value_adj numeric;
sequencevar text := gridtable ||'_ogc_fid_seq';
BEGIN
-- Check to see if grid table already exists
SELECT COUNT(*) INTO tbl_cnt FROM information_schema.tables WHERE table_schema = schemaname AND table_name = gridtable;
-- If grid table already exists, drop table
IF (tbl_cnt > 0) THEN
EXECUTE 'DROP TABLE ' || gridtable;
ELSE
END IF;
-- Create grid table
EXECUTE 'CREATE TABLE ' || gridtable ||'
(
ogc_fid serial,
x_count integer,
x numeric(8,1),
y_count integer,
y numeric(8,1),
geometry geometry(Polygon,3857)
)
WITH (
OIDS=FALSE
);
CREATE INDEX ' || gridtable || '_geom
ON ' || gridtable ||'
USING gist
(geometry);';
-- load extents of the bounding box and values for the initial y parameters
EXECUTE 'SELECT floor(st_xmin(ST_transform(geometry,3857))) FROM '|| boundingbox INTO XMIN;
EXECUTE 'SELECT ceiling(st_xmax(ST_transform(geometry,3857))) FROM '|| boundingbox INTO XMAX;
EXECUTE 'SELECT floor(st_ymin(ST_transform(geometry,3857))) FROM '|| boundingbox INTO YMIN;
EXECUTE 'SELECT ceiling(st_ymax(ST_transform(geometry,3857))) FROM '|| boundingbox INTO YMAX;
y_count = 1;
y_value = YMIN;
LOOP
-- for each y value, reset x to XMIN and subloop through the x values
x_count = 1;
x_value = XMIN;
LOOP
-- for every even numbered x_count, offset y by half the height of the hexagon
y_offset = (ceiling((x_count+1)/2::numeric) - floor((x_count+1)/2::numeric))*round(halfwidth*SQRT(3)::numeric,1)/2;
-- add the offset to the y_value
y_value_adj = y_value + y_offset;
EXECUTE 'INSERT INTO ' || gridtable || ' VALUES(NEXTVAL('|| quote_literal(sequencevar) ||'), '|| x_count ||', '|| x_value ||', ' || y_count ||', ' || y_value_adj ||', NULL)';
x_count = x_count + 1;
x_value = x_value + round(halfwidth*3/2::numeric,1);
EXIT WHEN x_value > XMAX;
END LOOP;
-- after exiting the subloop, increment the y count and y value
y_count = y_count + 1;
y_value = y_value + round(halfwidth*SQRT(3)::numeric,1);
EXIT WHEN y_value > YMAX;
END LOOP;
-- With the table now populated with x y points, the last step is to create a hexagon geometry for each.
EXECUTE 'UPDATE '|| gridtable ||' SET geometry = ST_transform(ST_Polygon(ST_makeline(ARRAY[st_makepoint(x-(0.5*'|| halfwidth ||'), y+(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x-(1*'|| halfwidth ||'), y), st_makepoint(x-(0.5*'|| halfwidth ||'), y-(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x+(0.5*'|| halfwidth ||'), y-(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x+(1*'|| halfwidth ||'), y), st_makepoint(x+(0.5*'|| halfwidth ||'), y+(SQRT(3)*0.5*'|| halfwidth ||')), st_makepoint(x-(0.5*'|| halfwidth ||'), y+(SQRT(3)*0.5*'|| halfwidth ||'))]),3857),3857)';
RETURN 'HEXAGON GRID CREATED';
END;
$BODY$
LANGUAGE plpgsql VOLATILE STRICT
COST 100;
-- USAGE:
-- SELECT makegrid_epsg3857('public',E'(SELECT ST_SetSRID(ST_MakePolygon(ST_GeomFromText(\'LINESTRING(-19.6875 26.1949,5.7129 26.1949,5.7129 44.2137,-19.6875 44.2137,-19.6875 26.1949)\')),4326) as geometry) a', 'hexamap', 10750)
-- REFERENCES:
-- http://www.dimensionaledge.com/main/postgis/how-to-create-hexagonal-grids-in-postgis/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment