pd.if.org Git - hexagon/blob - pgsql/hexfuncs.sql

   1 set role gis;
   2
   3 drop schema if exists hex cascade;
   4 create schema hex;
   5
   6 set search_path to hex,public,gis;
   7
   8 create or replace function hexgrid(
   9         xdim integer,
  10         ydim integer,
  11         angle numeric,
  12         coverage geometry,
  13         srid    integer,
  14         origin  geometry(POINT),
  15         tilesize numeric
  16 ) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
  17 declare
  18         rcover geometry;
  19         stride double precision;
  20         x integer;
  21         y integer;
  22         xmin double precision;
  23         ymin double precision;
  24         xoff double precision;
  25 begin
  26         -- the origin will be the center of the 1,1 hex.
  27         -- if null, then the origin will be set to fit in the
  28         -- covering box.
  29
  30         -- xdim, and ydim can be null, in which case they are
  31         -- calculated from the tilesize and the covering box
  32
  33         -- calculate ydim
  34
  35         -- precedence: ydim, xdim, tilesize
  36
  37         -- counter rotate geometry, fill, rotate grid
  38         rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
  39
  40         xmin := ST_XMin(rcover);
  41         ymin := ST_YMin(rcover);
  42
  43         stride := tilesize * 2.0 / sqrt(3.0);
  44
  45         --raise notice 'stride = %', stride;
  46         xdim := ((ST_XMax(rcover) -  ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
  47         ydim := ((ST_YMax(rcover) -  ST_YMin(rcover))/tilesize)::integer;
  48
  49         --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
  50
  51         -- TODO reorient these so right hand rule is inside
  52         hex := ST_Scale(ST_GeomFromText('POLYGON((
  53 .577350269189625764509148780502 0.0,
  54 .288675134594812882254574390251 0.5,
  55  -.288675134594812882254574390251 0.5,
  56  -.577350269189625764509148780502 0.0,
  57  -.288675134594812882254574390251 -0.5,
  58  .288675134594812882254574390251 -0.5,
  59  .577350269189625764509148780502 0.0))',ST_SRID(rcover)),tilesize,tilesize,1.0);
  60
  61         xoff := .577350269189625764509148780502 * 1.5;
  62
  63         for x in 0 .. xdim loop
  64                 for y in 0 .. ydim loop
  65                         -- TODO actually calculate cantor
  66                         return query select x, ydim-y, hex.cantor(ARRAY[x,ydim-y]),
  67                                 ST_Translate(hex,
  68                                 tilesize * (x * xoff) + xmin,
  69                                 --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
  70                                 tilesize * (y + 0.5 * (x%2) ) + ymin
  71                                 )
  72                                 --st_rotate(ST_Translate(hex,
  73                                 --tilesize * (x * xoff) + xmin,
  74                                 --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
  75                                 --), -3.1415926/6)
  76                         ;
  77                 end loop;
  78         end loop;
  79 END;
  80 $$
  81 language 'plpgsql'
  82 ;
  83
  84 -- create some hexagon ops
  85
  86 create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
  87 declare
  88         ijk integer[];
  89 begin
  90         if array_length(xy,1) = 3 then
  91                 return xy;
  92         end if;
  93
  94         ijk[1] := xy[1];
  95         ijk[2] := -xy[2];
  96
  97         if xy[1] < 0 then
  98                 ijk[2] := ijk[2] + (-xy[1]+1) / 2;
  99         else
 100                 ijk[2] := ijk[2] - xy[1]/2;
 101         end if;
 102
 103         ijk[3] := -ijk[1] - ijk[2];
 104
 105         return ijk;
 106 end
 107 $$ immutable strict;
 108
 109 create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
 110         select hex.ijk(ARRAY[x,y]);
 111 $$ immutable strict;
 112
 113 -- invert the cantor paring
 114 create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
 115 declare
 116         xy integer[];
 117         w integer;
 118         t integer;
 119         py integer;
 120 begin
 121         if hex >= 0 then
 122                 hex := hex-1;
 123                 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
 124                 t := (w * w + w)/2;
 125                 py := hex - t;
 126
 127                 xy = ARRAY[w-py,py];
 128         -- handle the negative cases
 129         else
 130                 hex := -hex;
 131                 hex := hex-1;
 132                 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
 133                 t := (w * w + w)/2;
 134                 py := hex - t;
 135
 136                 xy = ARRAY[w-py,py];
 137
 138                 if xy[1] % 2 = 1 then
 139                         xy[1] = -(xy[1]+1)/2;
 140                 else
 141                         xy[1] = xy[1]/2;
 142                 end if;
 143                 if xy[2] % 2 = 1 then
 144                         xy[2] = -(xy[2]+1)/2;
 145                 else
 146                         xy[2] = xy[2]/2;
 147                 end if;
 148         end if;
 149
 150         return xy;
 151 end;
 152 $$ immutable strict;
 153
 154 create function ijk(hex numeric) returns integer[] language 'sql' as $$
 155         select hex.ijk(hex.xy(hex));
 156 $$ immutable strict;
 157
 158 create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
 159 declare
 160         xy integer[];
 161 begin
 162         -- make sure this is really an ijk
 163         -- just no-op if it's already an
 164         -- xy, so this is sort of a coercion function
 165         if array_length(ijk,1) = 2 then
 166                 return ijk;
 167         end if;
 168         xy[1] := ijk[1];
 169         xy[2] := -ijk[2];
 170
 171         if ijk[1] < 0 then
 172                 xy[2] := xy[2] + (-ijk[1] + 1) / 2;
 173         else
 174                 xy[2] := xy[2] - ijk[1]/2;
 175         end if;
 176         return xy;
 177 end;
 178 $$;
 179
 180 CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
 181 RETURNS integer[] AS
 182 $$
 183 SELECT array_agg(result)
 184 FROM (SELECT tuple.val1 + tuple.val2 AS result
 185       FROM (SELECT UNNEST($1) AS val1
 186                    ,UNNEST($2) AS val2
 187                    ,generate_subscripts($1, 1) AS ix) tuple
 188       ORDER BY ix) inn;
 189 $$ LANGUAGE SQL IMMUTABLE STRICT;
 190
 191 create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
 192 declare
 193         adj integer[];
 194 begin
 195         if array_length(hex,1) = 2 then
 196                 hex = hex.ijk(hex);
 197         end if;
 198
 199         dir := dir % 6;
 200         if dir = 2 then
 201                 hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
 202         elsif dir = 1 then
 203                 hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
 204         elsif dir = 0 then
 205                 hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
 206         elsif dir = 5 then
 207                 hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
 208         elsif dir = 4 then
 209                 hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
 210         elsif dir = 3 then
 211                 hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
 212         else
 213                 -- error
 214         end if;
 215         return hex;
 216 end;
 217 $$ immutable strict;
 218
 219 create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
 220         select hex.adjhex(dir, hex.ijk(hex));
 221 $$ immutable strict;
 222
 223 -- TODO make this work for any, or inline it
 224 create function natcantor(hex integer[]) returns numeric language 'sql' as $$
 225         select cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 226 $$ immutable strict;
 227
 228 create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
 229 begin
 230         -- convert to xy
 231         if array_length(hex,1) = 3 then
 232                 hex = hex.xy(hex);
 233         end if;
 234
 235         if hex[1] < 0 or hex[2] < 0 then
 236                 hex[1] = hex[1] * 2;
 237                 hex[2] = hex[2] * 2;
 238                 if hex[1] < 0 then
 239                         hex[1] = -hex[1] - 1;
 240                 end if;
 241                 if hex[2] < 0 then
 242                         hex[2] = -hex[2] - 1;
 243                 end if;
 244                 return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 245         end if;
 246
 247         return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 248 end;
 249 $$ immutable strict;
 250
 251 create function adjacent(hex numeric) returns setof numeric language 'sql' as $$
 252 select hex.cantor(hex.adjhex(dir, hex)) from generate_series(0,5) as dir;
 253 $$ immutable strict;
 254
 255 create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
 256         select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
 257 $$ immutable strict;
 258
 259 create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
 260 declare
 261         c integer[]; -- cantor coordinates
 262         base integer[];
 263 begin
 264         if range <= 0 then
 265                 return;
 266         end if;
 267         if include then
 268                 return next hex;
 269         end if;
 270
 271         -- we could in theory eliminate the inner loop and
 272         -- do a return query
 273         base = hex.ijk(hex);
 274         for r in 1 .. range loop
 275                 c = base;
 276                 c[1] := c[1] + r;
 277                 c[3] := -c[1] - c[2]; -- keep the invariant
 278                 hex := hex.cantor(c);
 279
 280                 for h in 0 .. r * 6 - 1 loop
 281                         return next hex;
 282                         hex = hex.adjacent(h/r+2,hex);
 283                 end loop;
 284         end loop;
 285         return;
 286 end;
 287 $$ immutable strict;
 288
 289 create function within_array(hex numeric, range integer = 1, include boolean = false) returns numeric[] language 'sql' as $$
 290 select array_agg(c) from hex.within(hex, range, include) c;
 291 $$ immutable strict;
 292
 293 create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
 294 declare
 295         c integer[]; -- cantor coordinates
 296 begin
 297         if range = 0 then
 298                 return next hex;
 299         end if;
 300         if range <= 0 then
 301                 return;
 302         end if;
 303
 304         c = ijk(hex);
 305         c[1] := c[1] + range;
 306         c[3] := -c[1] - c[2]; -- keep the invariant
 307         hex := hex.cantor(c);
 308
 309         for q in 0 .. range * 6 - 1 loop
 310                 return next hex;
 311                 hex = hex.adjacent(q/range+2,hex);
 312         end loop;
 313         return;
 314 end;
 315 $$ immutable strict;
 316
 317 select
 318 xy(ARRAY[1,2]),
 319 ijk(xy(ARRAY[1,2])),
 320 xy(ijk(xy(ARRAY[1,2])))
 321 ;
 322
 323 select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
 324 select
 325 xy(adjhex(dir, ijk(ARRAY[1,2]))),
 326 cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))),
 327 xy(cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))))
 328 from generate_series(0,5) as dir;
 329
 330 select cantor(ARRAY[4,3]);
 331 --select * from within(cantor(ARRAY[4,4]), 2);
 332 select h, xy(h) from within(cantor(ARRAY[4,3]), 2) as h;