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 scale, 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 $$
24 stride := scale * 2.0 / sqrt(3.0);
26 x[1] := (2 + 3 * xc ) * stride;
27 x[4] := (-2 + 3 * xc) * stride;
28 x[2] := (1 + 3 * xc ) * stride;
29 x[6] := (1 + 3 * xc ) * stride;
30 x[3] := (-1 + 3 * xc ) * stride;
31 x[5] := (-1 + 3 * xc ) * stride;
33 -- hex coordinates are positive down
34 y[1] := (0 - 2 * yc + xc % 2) * yoff;
35 y[4] := (0 - 2 * yc + xc % 2) * yoff;
36 y[2] := (1 - 2 * yc + xc % 2) * yoff;
37 y[3] := (1 - 2 * yc + xc % 2) * yoff;
38 y[5] := (-1 - 2 * yc + xc % 2) * yoff;
39 y[6] := (-1 - 2 * yc + xc % 2) * yoff;
46 xoff := cosa * x[i] + sina * y[i];
47 y[i] := -sina * x[i] + cosa * y[i];
51 if xorigin != 0.0 or yorigin != 0.0 then
53 y[i] := y[i] + yorigin;
54 x[i] := x[i] + xorigin;
58 return format('POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s, %s %s))',
59 x[1], y[1], x[2], y[2], x[3], y[3], x[4], y[4], x[5], y[5], x[6], y[6]);
63 create or replace function hexcovergrid(
68 srid integer, -- will default to the srid of the coverage
69 origin geometry(POINT) -- this would be the 0,0 hex center, or perhaps 1,1
70 ) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
76 stride double precision;
81 xoff double precision;
83 -- the origin will be the center of the 1,1 hex.
84 -- if null, then the origin will be set to fit in the
87 -- counter rotate geometry, fill, rotate grid
88 rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
89 if srid is null and rcover is not null then
90 srid := st_srid(rcover);
93 stride := hexsize * 2.0 / sqrt(3.0);
94 -- the center of the 1,1 hex should be one stride right and one half tilesize
95 -- down from the upper right of the coverage. That is, the upper right of
96 -- (the bounding box of) the coverage should coincide with the "northwest"
98 xorigin = st_xmin(rcover) + stride;
99 yorigin = st_ymax(rcover) + hexsize / 2.0;
101 xdim := ((ST_XMax(rcover) - ST_XMin(rcover))/hexsize*2/sqrt(3.0))::integer;
102 ydim := ((ST_YMax(rcover) - ST_YMin(rcover))/hexsize)::integer;
104 --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
106 for x in 0 .. xdim loop
107 for y in 0 .. ydim loop
108 hex := st_geomfromtext(hexgeom(x, ydim-y, hexsize,angle,xorigin,yorigin),srid);
109 if not full and not st_intersects(hex, rcover) then
114 hex.cantor(ARRAY[x,ydim-y]),
124 create or replace function hexgrid(
130 origin geometry(POINT),
132 ) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
135 stride double precision;
138 xmin double precision;
139 ymin double precision;
140 xoff double precision;
142 -- the origin will be the center of the 1,1 hex.
143 -- if null, then the origin will be set to fit in the
146 -- xdim, and ydim can be null, in which case they are
147 -- calculated from the tilesize and the covering box
151 -- precedence: ydim, xdim, tilesize
153 -- counter rotate geometry, fill, rotate grid
154 rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
155 if srid is null and rcover is not null then
156 srid := st_srid(rcover);
159 xmin := ST_XMin(rcover);
160 ymin := ST_YMin(rcover);
162 stride := tilesize * 2.0 / sqrt(3.0);
164 xdim := ((ST_XMax(rcover) - ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
165 ydim := ((ST_YMax(rcover) - ST_YMin(rcover))/tilesize)::integer;
167 --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
169 for x in 0 .. xdim loop
170 for y in 0 .. ydim loop
173 hex.cantor(ARRAY[x,ydim-y]),
174 st_geomfromtext(hexgeom(x, ydim-y, tilesize,angle),srid)
183 -- create some hexagon ops
185 create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
189 if array_length(xy,1) = 3 then
197 ijk[2] := ijk[2] + (-xy[1]+1) / 2;
199 ijk[2] := ijk[2] - xy[1]/2;
202 ijk[3] := -ijk[1] - ijk[2];
208 create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
209 select hex.ijk(ARRAY[x,y]);
212 -- invert the cantor paring
213 create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
222 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
227 -- handle the negative cases
231 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
237 if xy[1] % 2 = 1 then
238 xy[1] = -(xy[1]+1)/2;
242 if xy[2] % 2 = 1 then
243 xy[2] = -(xy[2]+1)/2;
253 create function ijk(hex numeric) returns integer[] language 'sql' as $$
254 select hex.ijk(hex.xy(hex));
257 create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
261 -- make sure this is really an ijk
262 -- just no-op if it's already an
263 -- xy, so this is sort of a coercion function
264 if array_length(ijk,1) = 2 then
271 xy[2] := xy[2] + (-ijk[1] + 1) / 2;
273 xy[2] := xy[2] - ijk[1]/2;
279 CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
282 SELECT array_agg(result)
283 FROM (SELECT tuple.val1 + tuple.val2 AS result
284 FROM (SELECT UNNEST($1) AS val1
286 ,generate_subscripts($1, 1) AS ix) tuple
288 $$ LANGUAGE SQL IMMUTABLE STRICT;
290 create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
294 if array_length(hex,1) = 2 then
300 hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
302 hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
304 hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
306 hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
308 hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
310 hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
318 create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
319 select hex.adjhex(dir, hex.ijk(hex));
322 -- TODO make this work for any, or inline it
323 create function natcantor(hex integer[]) returns numeric language 'sql' as $$
324 select cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
327 create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
330 if array_length(hex,1) = 3 then
334 if hex[1] < 0 or hex[2] < 0 then
338 hex[1] = -hex[1] - 1;
341 hex[2] = -hex[2] - 1;
343 return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
346 return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
350 create function adjacent(hex numeric) returns setof numeric language 'sql' as $$
351 select hex.cantor(hex.adjhex(dir, hex)) from generate_series(0,5) as dir;
354 create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
355 select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
358 create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
360 c integer[]; -- cantor coordinates
370 -- we could in theory eliminate the inner loop and
373 for r in 1 .. range loop
376 c[3] := -c[1] - c[2]; -- keep the invariant
377 hex := hex.cantor(c);
379 for h in 0 .. r * 6 - 1 loop
381 hex = hex.adjacent(h/r+2,hex);
388 create function within_array(hex numeric, range integer = 1, include boolean = false) returns numeric[] language 'sql' as $$
389 select array_agg(c) from hex.within(hex, range, include) c;
392 create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
394 c integer[]; -- cantor coordinates
404 c[1] := c[1] + range;
405 c[3] := -c[1] - c[2]; -- keep the invariant
406 hex := hex.cantor(c);
408 for q in 0 .. range * 6 - 1 loop
410 hex = hex.adjacent(q/range+2,hex);
419 xy(ijk(xy(ARRAY[1,2])))
422 select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
424 xy(adjhex(dir, ijk(ARRAY[1,2]))),
425 cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))),
426 xy(cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))))
427 from generate_series(0,5) as dir;
429 select cantor(ARRAY[4,3]);
430 --select * from within(cantor(ARRAY[4,4]), 2);
431 select h, xy(h) from within(cantor(ARRAY[4,3]), 2) as h;