+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'
+;
+