From 962564c7f79681915217cc92329f82a9a36442d3 Mon Sep 17 00:00:00 2001 From: Nathan Wagner Date: Tue, 17 Mar 2015 19:07:07 +0000 Subject: [PATCH] rewrite parts of postgres implementation --- pgsql/hexfuncs.sql | 128 +++++++++++++++++++++++++++++++++------------ 1 file changed, 96 insertions(+), 32 deletions(-) diff --git a/pgsql/hexfuncs.sql b/pgsql/hexfuncs.sql index b080539..ee8f349 100644 --- 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; -- 2.40.0