+-- return a postgis text polygon
+-- the 0,0 hex center is at the origin
+-- we scale, rotate, then translate
+create or replace function hexgeom(
+ xc integer, yc integer, scale float = 1.0, angle float = 0.0, xorigin float = 0.0, yorigin float = 0.0
+) returns text language 'plpgsql' as $$
+declare
+ x float[];
+ y float[];
+ stride float;
+ xoff float;
+ yoff float;
+ i integer;
+ cosa float;
+ sina float;
+begin
+ stride := scale * 2.0 / sqrt(3.0);
+ yoff = scale / 2.0;
+ x[1] := (2 + 3 * xc ) * stride;
+ x[4] := (-2 + 3 * xc) * stride;
+ x[2] := (1 + 3 * xc ) * stride;
+ x[6] := (1 + 3 * xc ) * stride;
+ x[3] := (-1 + 3 * xc ) * stride;
+ x[5] := (-1 + 3 * xc ) * stride;
+
+ -- hex coordinates are positive down
+ y[1] := (0 - 2 * yc + xc % 2) * yoff;
+ y[4] := (0 - 2 * yc + xc % 2) * yoff;
+ y[2] := (1 - 2 * yc + xc % 2) * yoff;
+ y[3] := (1 - 2 * yc + xc % 2) * yoff;
+ y[5] := (-1 - 2 * yc + xc % 2) * yoff;
+ y[6] := (-1 - 2 * yc + xc % 2) * yoff;
+
+ if angle != 0.0 then
+ cosa := cos(angle);
+ sina := sin(angle);
+
+ for i in 1..6 loop
+ xoff := cosa * x[i] + sina * y[i];
+ y[i] := -sina * x[i] + cosa * y[i];
+ x[i] := xoff;
+ end loop;
+ end if;
+ if xorigin != 0.0 or yorigin != 0.0 then
+ for i in 1..6 loop
+ y[i] := y[i] + yorigin;
+ x[i] := x[i] + xorigin;
+ end loop;
+ end if;
+
+ return format('POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s, %s %s))',
+ x[1], y[1], x[2], y[2], x[3], y[3], x[4], y[4], x[5], y[5], x[6], y[6]);
+end;
+$$ immutable strict;
+
+create or replace function hexcovergrid(
+ coverage geometry,
+ hexsize numeric
+ angle numeric = 0.0,
+ full boolean = true,
+ srid integer, -- will default to the srid of the coverage
+ origin geometry(POINT) -- this would be the 0,0 hex center, or perhaps 1,1
+) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
+declare
+ xdim integer,
+ ydim integer,
+ rcover geometry;
+ hex geometry;
+ stride double precision;
+ xorigin float;
+ yorigin float;
+ x integer;
+ y integer;
+ xoff double precision;
+begin
+ -- the origin will be the center of the 1,1 hex.
+ -- if null, then the origin will be set to fit in the
+ -- covering box.
+
+ -- counter rotate geometry, fill, rotate grid
+ rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
+ if srid is null and rcover is not null then
+ srid := st_srid(rcover);
+ end if;
+
+ stride := hexsize * 2.0 / sqrt(3.0);
+ -- the center of the 1,1 hex should be one stride right and one half tilesize
+ -- down from the upper right of the coverage. That is, the upper right of
+ -- (the bounding box of) the coverage should coincide with the "northwest"
+ -- vertex of hex 1,1
+ xorigin = st_xmin(rcover) + stride;
+ yorigin = st_ymax(rcover) + hexsize / 2.0;
+
+ xdim := ((ST_XMax(rcover) - ST_XMin(rcover))/hexsize*2/sqrt(3.0))::integer;
+ ydim := ((ST_YMax(rcover) - ST_YMin(rcover))/hexsize)::integer;
+
+ --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
+
+ for x in 0 .. xdim loop
+ for y in 0 .. ydim loop
+ hex := st_geomfromtext(hexgeom(x, ydim-y, hexsize,angle,xorigin,yorigin),srid);
+ if not full and not st_intersects(hex, rcover) then
+ next;
+ end if;
+ return query select
+ x, ydim-y,
+ hex.cantor(ARRAY[x,ydim-y]),
+ hex
+ ;
+ end loop;
+ end loop;
+END;
+$$
+language 'plpgsql'
+;
+