rewrite parts of postgres implementation

diff --git a/pgsql/hexfuncs.sql b/pgsql/hexfuncs.sql
index b080539cac71a03a54f1bece2ed4c417b971bf89..ee8f34959dfc064f43a486523130621b71ab0fa3 100644 (file)
--- a/pgsql/hexfuncs.sql
+++ b/pgsql/hexfuncs.sql
@@ -7,7 +7,7 @@ set search_path to hex,public,gis;
 
 -- return a postgis text polygon
 -- the 0,0 hex center is at the origin
--- we rotate, then translate
+-- 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 $$
@@ -17,9 +17,11 @@ declare
        stride float;
        xoff float;
        yoff float;
+       i integer;
+       cosa float;
+       sina float;
 begin
        stride := scale * 2.0 / sqrt(3.0);
-       xoff := stride;
        yoff = scale / 2.0;
        x[1] := (2 + 3 * xc ) * stride;
        x[4] := (-2 + 3 * xc) * stride;
@@ -28,18 +30,97 @@ begin
        x[3] := (-1 + 3 * xc ) * stride;
        x[5] := (-1 + 3 * xc ) * stride;
 
-       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;
+       -- 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'
+;
+
 create or replace function hexgrid(
        xdim integer,
        ydim integer,
@@ -71,43 +152,26 @@ begin
 
        -- 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;
        
        xmin := ST_XMin(rcover);
        ymin := ST_YMin(rcover);
        
        stride := tilesize * 2.0 / sqrt(3.0);
 
-       --raise notice 'stride = %', stride;
        xdim := ((ST_XMax(rcover) -  ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
        ydim := ((ST_YMax(rcover) -  ST_YMin(rcover))/tilesize)::integer;
 
        --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
        
-        -- TODO reorient these so right hand rule is inside
-        hex := ST_Scale(ST_GeomFromText('POLYGON((
-.577350269189625764509148780502 0.0,
-.288675134594812882254574390251 0.5,
- -.288675134594812882254574390251 0.5,
- -.577350269189625764509148780502 0.0,
- -.288675134594812882254574390251 -0.5,
- .288675134594812882254574390251 -0.5,
- .577350269189625764509148780502 0.0))',ST_SRID(rcover)),tilesize,tilesize,1.0);
-
-        xoff := .577350269189625764509148780502 * 1.5;
-
         for x in 0 .. xdim loop
                 for y in 0 .. ydim loop
-                       -- TODO actually calculate cantor
-                       return query select x, ydim-y, hex.cantor(ARRAY[x,ydim-y]), 
-                               ST_Translate(hex,
-                               tilesize * (x * xoff) + xmin,
-                               --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
-                               tilesize * (y + 0.5 * (x%2) ) + ymin
-                               )
-                               --st_rotate(ST_Translate(hex,
-                               --tilesize * (x * xoff) + xmin,
-                               --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
-                               --), -3.1415926/6)
+                       return query select
+                               x, ydim-y,
+                               hex.cantor(ARRAY[x,ydim-y]), 
+                               st_geomfromtext(hexgeom(x, ydim-y, tilesize,angle),srid)
                        ;
                 end loop;
         end loop;