3 drop schema if exists hex cascade;
6 set search_path to hex,public,gis;
8 -- return a postgis text polygon
9 -- the 0,0 hex center is at the origin
10 -- we rotate, then translate
11 create or replace function hexgeom(
12 xc integer, yc integer, scale float = 1.0, angle float = 0.0, xorigin float = 0.0, yorigin float = 0.0
13 ) returns text language 'plpgsql' as $$
21 stride := scale * 2.0 / sqrt(3.0);
24 x[1] := (2 + 3 * xc ) * stride;
25 x[4] := (-2 + 3 * xc) * stride;
26 x[2] := (1 + 3 * xc ) * stride;
27 x[6] := (1 + 3 * xc ) * stride;
28 x[3] := (-1 + 3 * xc ) * stride;
29 x[5] := (-1 + 3 * xc ) * stride;
31 y[1] := (0 + 2 * yc + xc % 2) * yoff;
32 y[4] := (0 + 2 * yc + xc % 2) * yoff;
33 y[2] := (1 + 2 * yc + xc % 2) * yoff;
34 y[3] := (1 + 2 * yc + xc % 2) * yoff;
35 y[5] := (-1 + 2 * yc + xc % 2) * yoff;
36 y[6] := (-1 + 2 * yc + xc % 2) * yoff;
38 return format('POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s, %s %s))',
39 x[1], y[1], x[2], y[2], x[3], y[3], x[4], y[4], x[5], y[5], x[6], y[6]);
43 create or replace function hexgrid(
49 origin geometry(POINT),
51 ) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
54 stride double precision;
57 xmin double precision;
58 ymin double precision;
59 xoff double precision;
61 -- the origin will be the center of the 1,1 hex.
62 -- if null, then the origin will be set to fit in the
65 -- xdim, and ydim can be null, in which case they are
66 -- calculated from the tilesize and the covering box
70 -- precedence: ydim, xdim, tilesize
72 -- counter rotate geometry, fill, rotate grid
73 rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
75 xmin := ST_XMin(rcover);
76 ymin := ST_YMin(rcover);
78 stride := tilesize * 2.0 / sqrt(3.0);
80 --raise notice 'stride = %', stride;
81 xdim := ((ST_XMax(rcover) - ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
82 ydim := ((ST_YMax(rcover) - ST_YMin(rcover))/tilesize)::integer;
84 --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
86 -- TODO reorient these so right hand rule is inside
87 hex := ST_Scale(ST_GeomFromText('POLYGON((
88 .577350269189625764509148780502 0.0,
89 .288675134594812882254574390251 0.5,
90 -.288675134594812882254574390251 0.5,
91 -.577350269189625764509148780502 0.0,
92 -.288675134594812882254574390251 -0.5,
93 .288675134594812882254574390251 -0.5,
94 .577350269189625764509148780502 0.0))',ST_SRID(rcover)),tilesize,tilesize,1.0);
96 xoff := .577350269189625764509148780502 * 1.5;
98 for x in 0 .. xdim loop
99 for y in 0 .. ydim loop
100 -- TODO actually calculate cantor
101 return query select x, ydim-y, hex.cantor(ARRAY[x,ydim-y]),
103 tilesize * (x * xoff) + xmin,
104 --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
105 tilesize * (y + 0.5 * (x%2) ) + ymin
107 --st_rotate(ST_Translate(hex,
108 --tilesize * (x * xoff) + xmin,
109 --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
119 -- create some hexagon ops
121 create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
125 if array_length(xy,1) = 3 then
133 ijk[2] := ijk[2] + (-xy[1]+1) / 2;
135 ijk[2] := ijk[2] - xy[1]/2;
138 ijk[3] := -ijk[1] - ijk[2];
144 create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
145 select hex.ijk(ARRAY[x,y]);
148 -- invert the cantor paring
149 create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
158 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
163 -- handle the negative cases
167 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
173 if xy[1] % 2 = 1 then
174 xy[1] = -(xy[1]+1)/2;
178 if xy[2] % 2 = 1 then
179 xy[2] = -(xy[2]+1)/2;
189 create function ijk(hex numeric) returns integer[] language 'sql' as $$
190 select hex.ijk(hex.xy(hex));
193 create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
197 -- make sure this is really an ijk
198 -- just no-op if it's already an
199 -- xy, so this is sort of a coercion function
200 if array_length(ijk,1) = 2 then
207 xy[2] := xy[2] + (-ijk[1] + 1) / 2;
209 xy[2] := xy[2] - ijk[1]/2;
215 CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
218 SELECT array_agg(result)
219 FROM (SELECT tuple.val1 + tuple.val2 AS result
220 FROM (SELECT UNNEST($1) AS val1
222 ,generate_subscripts($1, 1) AS ix) tuple
224 $$ LANGUAGE SQL IMMUTABLE STRICT;
226 create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
230 if array_length(hex,1) = 2 then
236 hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
238 hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
240 hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
242 hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
244 hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
246 hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
254 create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
255 select hex.adjhex(dir, hex.ijk(hex));
258 -- TODO make this work for any, or inline it
259 create function natcantor(hex integer[]) returns numeric language 'sql' as $$
260 select cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
263 create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
266 if array_length(hex,1) = 3 then
270 if hex[1] < 0 or hex[2] < 0 then
274 hex[1] = -hex[1] - 1;
277 hex[2] = -hex[2] - 1;
279 return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
282 return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
286 create function adjacent(hex numeric) returns setof numeric language 'sql' as $$
287 select hex.cantor(hex.adjhex(dir, hex)) from generate_series(0,5) as dir;
290 create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
291 select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
294 create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
296 c integer[]; -- cantor coordinates
306 -- we could in theory eliminate the inner loop and
309 for r in 1 .. range loop
312 c[3] := -c[1] - c[2]; -- keep the invariant
313 hex := hex.cantor(c);
315 for h in 0 .. r * 6 - 1 loop
317 hex = hex.adjacent(h/r+2,hex);
324 create function within_array(hex numeric, range integer = 1, include boolean = false) returns numeric[] language 'sql' as $$
325 select array_agg(c) from hex.within(hex, range, include) c;
328 create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
330 c integer[]; -- cantor coordinates
340 c[1] := c[1] + range;
341 c[3] := -c[1] - c[2]; -- keep the invariant
342 hex := hex.cantor(c);
344 for q in 0 .. range * 6 - 1 loop
346 hex = hex.adjacent(q/range+2,hex);
355 xy(ijk(xy(ARRAY[1,2])))
358 select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
360 xy(adjhex(dir, ijk(ARRAY[1,2]))),
361 cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))),
362 xy(cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))))
363 from generate_series(0,5) as dir;
365 select cantor(ARRAY[4,3]);
366 --select * from within(cantor(ARRAY[4,4]), 2);
367 select h, xy(h) from within(cantor(ARRAY[4,3]), 2) as h;