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 -- 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 $$
  14 declare
  15         x float[];
  16         y float[];
  17         stride float;
  18         xoff float;
  19         yoff float;
  20         i integer;
  21         cosa float;
  22         sina float;
  23 begin
  24         stride := scale * 2.0 / sqrt(3.0);
  25         yoff = scale / 2.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;
  32
  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;
  40
  41         if angle != 0.0 then
  42                 cosa := cos(angle);
  43                 sina := sin(angle);
  44
  45                 for i in 1..6 loop
  46                         xoff := cosa * x[i] + sina * y[i];
  47                         y[i] := -sina * x[i] + cosa * y[i];
  48                         x[i] := xoff;
  49                 end loop;
  50         end if;
  51         if xorigin != 0.0 or yorigin != 0.0 then
  52                 for i in 1..6 loop
  53                         y[i] := y[i] + yorigin;
  54                         x[i] := x[i] + xorigin;
  55                 end loop;
  56         end if;
  57
  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]);
  60 end;
  61 $$ immutable strict;
  62
  63 create or replace function hexcovergrid(
  64         coverage geometry,
  65         hexsize numeric
  66         angle numeric = 0.0,
  67         full boolean = true,
  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 $$
  71 declare
  72         xdim integer,
  73         ydim integer,
  74         rcover geometry;
  75         hex geometry;
  76         stride double precision;
  77         xorigin float;
  78         yorigin float;
  79         x integer;
  80         y integer;
  81         xoff double precision;
  82 begin
  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
  85         -- covering box.
  86
  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);
  91         end if;
  92
  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"
  97         -- vertex of hex 1,1
  98         xorigin = st_xmin(rcover) + stride;
  99         yorigin = st_ymax(rcover) + hexsize / 2.0;
 100
 101         xdim := ((ST_XMax(rcover) -  ST_XMin(rcover))/hexsize*2/sqrt(3.0))::integer;
 102         ydim := ((ST_YMax(rcover) -  ST_YMin(rcover))/hexsize)::integer;
 103
 104         --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
 105
 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
 110                                 next;
 111                         end if;
 112                         return query select
 113                                 x, ydim-y,
 114                                 hex.cantor(ARRAY[x,ydim-y]),
 115                                 hex
 116                         ;
 117                 end loop;
 118         end loop;
 119 END;
 120 $$
 121 language 'plpgsql'
 122 ;
 123
 124 create or replace function hexgrid(
 125         xdim integer,
 126         ydim integer,
 127         angle numeric,
 128         coverage geometry,
 129         srid    integer,
 130         origin  geometry(POINT),
 131         tilesize numeric
 132 ) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
 133 declare
 134         rcover geometry;
 135         stride double precision;
 136         x integer;
 137         y integer;
 138         xmin double precision;
 139         ymin double precision;
 140         xoff double precision;
 141 begin
 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
 144         -- covering box.
 145
 146         -- xdim, and ydim can be null, in which case they are
 147         -- calculated from the tilesize and the covering box
 148
 149         -- calculate ydim
 150
 151         -- precedence: ydim, xdim, tilesize
 152
 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);
 157         end if;
 158
 159         xmin := ST_XMin(rcover);
 160         ymin := ST_YMin(rcover);
 161
 162         stride := tilesize * 2.0 / sqrt(3.0);
 163
 164         xdim := ((ST_XMax(rcover) -  ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
 165         ydim := ((ST_YMax(rcover) -  ST_YMin(rcover))/tilesize)::integer;
 166
 167         --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
 168
 169         for x in 0 .. xdim loop
 170                 for y in 0 .. ydim loop
 171                         return query select
 172                                 x, ydim-y,
 173                                 hex.cantor(ARRAY[x,ydim-y]),
 174                                 st_geomfromtext(hexgeom(x, ydim-y, tilesize,angle),srid)
 175                         ;
 176                 end loop;
 177         end loop;
 178 END;
 179 $$
 180 language 'plpgsql'
 181 ;
 182
 183 -- create some hexagon ops
 184
 185 create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
 186 declare
 187         ijk integer[];
 188 begin
 189         if array_length(xy,1) = 3 then
 190                 return xy;
 191         end if;
 192
 193         ijk[1] := xy[1];
 194         ijk[2] := -xy[2];
 195
 196         if xy[1] < 0 then
 197                 ijk[2] := ijk[2] + (-xy[1]+1) / 2;
 198         else
 199                 ijk[2] := ijk[2] - xy[1]/2;
 200         end if;
 201
 202         ijk[3] := -ijk[1] - ijk[2];
 203
 204         return ijk;
 205 end
 206 $$ immutable strict;
 207
 208 create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
 209         select hex.ijk(ARRAY[x,y]);
 210 $$ immutable strict;
 211
 212 -- invert the cantor paring
 213 create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
 214 declare
 215         xy integer[];
 216         w integer;
 217         t integer;
 218         py integer;
 219 begin
 220         if hex >= 0 then
 221                 hex := hex-1;
 222                 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
 223                 t := (w * w + w)/2;
 224                 py := hex - t;
 225
 226                 xy = ARRAY[w-py,py];
 227         -- handle the negative cases
 228         else
 229                 hex := -hex;
 230                 hex := hex-1;
 231                 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
 232                 t := (w * w + w)/2;
 233                 py := hex - t;
 234
 235                 xy = ARRAY[w-py,py];
 236
 237                 if xy[1] % 2 = 1 then
 238                         xy[1] = -(xy[1]+1)/2;
 239                 else
 240                         xy[1] = xy[1]/2;
 241                 end if;
 242                 if xy[2] % 2 = 1 then
 243                         xy[2] = -(xy[2]+1)/2;
 244                 else
 245                         xy[2] = xy[2]/2;
 246                 end if;
 247         end if;
 248
 249         return xy;
 250 end;
 251 $$ immutable strict;
 252
 253 create function ijk(hex numeric) returns integer[] language 'sql' as $$
 254         select hex.ijk(hex.xy(hex));
 255 $$ immutable strict;
 256
 257 create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
 258 declare
 259         xy integer[];
 260 begin
 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
 265                 return ijk;
 266         end if;
 267         xy[1] := ijk[1];
 268         xy[2] := -ijk[2];
 269
 270         if ijk[1] < 0 then
 271                 xy[2] := xy[2] + (-ijk[1] + 1) / 2;
 272         else
 273                 xy[2] := xy[2] - ijk[1]/2;
 274         end if;
 275         return xy;
 276 end;
 277 $$;
 278
 279 CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
 280 RETURNS integer[] AS
 281 $$
 282 SELECT array_agg(result)
 283 FROM (SELECT tuple.val1 + tuple.val2 AS result
 284       FROM (SELECT UNNEST($1) AS val1
 285                    ,UNNEST($2) AS val2
 286                    ,generate_subscripts($1, 1) AS ix) tuple
 287       ORDER BY ix) inn;
 288 $$ LANGUAGE SQL IMMUTABLE STRICT;
 289
 290 create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
 291 declare
 292         adj integer[];
 293 begin
 294         if array_length(hex,1) = 2 then
 295                 hex = hex.ijk(hex);
 296         end if;
 297
 298         dir := dir % 6;
 299         if dir = 2 then
 300                 hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
 301         elsif dir = 1 then
 302                 hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
 303         elsif dir = 0 then
 304                 hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
 305         elsif dir = 5 then
 306                 hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
 307         elsif dir = 4 then
 308                 hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
 309         elsif dir = 3 then
 310                 hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
 311         else
 312                 -- error
 313         end if;
 314         return hex;
 315 end;
 316 $$ immutable strict;
 317
 318 create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
 319         select hex.adjhex(dir, hex.ijk(hex));
 320 $$ immutable strict;
 321
 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);
 325 $$ immutable strict;
 326
 327 create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
 328 begin
 329         -- convert to xy
 330         if array_length(hex,1) = 3 then
 331                 hex = hex.xy(hex);
 332         end if;
 333
 334         if hex[1] < 0 or hex[2] < 0 then
 335                 hex[1] = hex[1] * 2;
 336                 hex[2] = hex[2] * 2;
 337                 if hex[1] < 0 then
 338                         hex[1] = -hex[1] - 1;
 339                 end if;
 340                 if hex[2] < 0 then
 341                         hex[2] = -hex[2] - 1;
 342                 end if;
 343                 return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 344         end if;
 345
 346         return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 347 end;
 348 $$ immutable strict;
 349
 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;
 352 $$ immutable strict;
 353
 354 create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
 355         select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
 356 $$ immutable strict;
 357
 358 create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
 359 declare
 360         c integer[]; -- cantor coordinates
 361         base integer[];
 362 begin
 363         if range <= 0 then
 364                 return;
 365         end if;
 366         if include then
 367                 return next hex;
 368         end if;
 369
 370         -- we could in theory eliminate the inner loop and
 371         -- do a return query
 372         base = hex.ijk(hex);
 373         for r in 1 .. range loop
 374                 c = base;
 375                 c[1] := c[1] + r;
 376                 c[3] := -c[1] - c[2]; -- keep the invariant
 377                 hex := hex.cantor(c);
 378
 379                 for h in 0 .. r * 6 - 1 loop
 380                         return next hex;
 381                         hex = hex.adjacent(h/r+2,hex);
 382                 end loop;
 383         end loop;
 384         return;
 385 end;
 386 $$ immutable strict;
 387
 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;
 390 $$ immutable strict;
 391
 392 create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
 393 declare
 394         c integer[]; -- cantor coordinates
 395 begin
 396         if range = 0 then
 397                 return next hex;
 398         end if;
 399         if range <= 0 then
 400                 return;
 401         end if;
 402
 403         c = ijk(hex);
 404         c[1] := c[1] + range;
 405         c[3] := -c[1] - c[2]; -- keep the invariant
 406         hex := hex.cantor(c);
 407
 408         for q in 0 .. range * 6 - 1 loop
 409                 return next hex;
 410                 hex = hex.adjacent(q/range+2,hex);
 411         end loop;
 412         return;
 413 end;
 414 $$ immutable strict;
 415
 416 select
 417 xy(ARRAY[1,2]),
 418 ijk(xy(ARRAY[1,2])),
 419 xy(ijk(xy(ARRAY[1,2])))
 420 ;
 421
 422 select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
 423 select
 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;
 428
 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;