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 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 begin
  21         stride := scale * 2.0 / sqrt(3.0);
  22         xoff := stride;
  23         yoff = scale / 2.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;
  30
  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;
  37
  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]);
  40 end;
  41 $$ immutable strict;
  42
  43 create or replace function hexgrid(
  44         xdim integer,
  45         ydim integer,
  46         angle numeric,
  47         coverage geometry,
  48         srid    integer,
  49         origin  geometry(POINT),
  50         tilesize numeric
  51 ) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
  52 declare
  53         rcover geometry;
  54         stride double precision;
  55         x integer;
  56         y integer;
  57         xmin double precision;
  58         ymin double precision;
  59         xoff double precision;
  60 begin
  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
  63         -- covering box.
  64
  65         -- xdim, and ydim can be null, in which case they are
  66         -- calculated from the tilesize and the covering box
  67
  68         -- calculate ydim
  69
  70         -- precedence: ydim, xdim, tilesize
  71
  72         -- counter rotate geometry, fill, rotate grid
  73         rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
  74
  75         xmin := ST_XMin(rcover);
  76         ymin := ST_YMin(rcover);
  77
  78         stride := tilesize * 2.0 / sqrt(3.0);
  79
  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;
  83
  84         --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
  85
  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);
  95
  96         xoff := .577350269189625764509148780502 * 1.5;
  97
  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]),
 102                                 ST_Translate(hex,
 103                                 tilesize * (x * xoff) + xmin,
 104                                 --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
 105                                 tilesize * (y + 0.5 * (x%2) ) + ymin
 106                                 )
 107                                 --st_rotate(ST_Translate(hex,
 108                                 --tilesize * (x * xoff) + xmin,
 109                                 --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
 110                                 --), -3.1415926/6)
 111                         ;
 112                 end loop;
 113         end loop;
 114 END;
 115 $$
 116 language 'plpgsql'
 117 ;
 118
 119 -- create some hexagon ops
 120
 121 create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
 122 declare
 123         ijk integer[];
 124 begin
 125         if array_length(xy,1) = 3 then
 126                 return xy;
 127         end if;
 128
 129         ijk[1] := xy[1];
 130         ijk[2] := -xy[2];
 131
 132         if xy[1] < 0 then
 133                 ijk[2] := ijk[2] + (-xy[1]+1) / 2;
 134         else
 135                 ijk[2] := ijk[2] - xy[1]/2;
 136         end if;
 137
 138         ijk[3] := -ijk[1] - ijk[2];
 139
 140         return ijk;
 141 end
 142 $$ immutable strict;
 143
 144 create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
 145         select hex.ijk(ARRAY[x,y]);
 146 $$ immutable strict;
 147
 148 -- invert the cantor paring
 149 create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
 150 declare
 151         xy integer[];
 152         w integer;
 153         t integer;
 154         py integer;
 155 begin
 156         if hex >= 0 then
 157                 hex := hex-1;
 158                 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
 159                 t := (w * w + w)/2;
 160                 py := hex - t;
 161
 162                 xy = ARRAY[w-py,py];
 163         -- handle the negative cases
 164         else
 165                 hex := -hex;
 166                 hex := hex-1;
 167                 w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
 168                 t := (w * w + w)/2;
 169                 py := hex - t;
 170
 171                 xy = ARRAY[w-py,py];
 172
 173                 if xy[1] % 2 = 1 then
 174                         xy[1] = -(xy[1]+1)/2;
 175                 else
 176                         xy[1] = xy[1]/2;
 177                 end if;
 178                 if xy[2] % 2 = 1 then
 179                         xy[2] = -(xy[2]+1)/2;
 180                 else
 181                         xy[2] = xy[2]/2;
 182                 end if;
 183         end if;
 184
 185         return xy;
 186 end;
 187 $$ immutable strict;
 188
 189 create function ijk(hex numeric) returns integer[] language 'sql' as $$
 190         select hex.ijk(hex.xy(hex));
 191 $$ immutable strict;
 192
 193 create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
 194 declare
 195         xy integer[];
 196 begin
 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
 201                 return ijk;
 202         end if;
 203         xy[1] := ijk[1];
 204         xy[2] := -ijk[2];
 205
 206         if ijk[1] < 0 then
 207                 xy[2] := xy[2] + (-ijk[1] + 1) / 2;
 208         else
 209                 xy[2] := xy[2] - ijk[1]/2;
 210         end if;
 211         return xy;
 212 end;
 213 $$;
 214
 215 CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
 216 RETURNS integer[] AS
 217 $$
 218 SELECT array_agg(result)
 219 FROM (SELECT tuple.val1 + tuple.val2 AS result
 220       FROM (SELECT UNNEST($1) AS val1
 221                    ,UNNEST($2) AS val2
 222                    ,generate_subscripts($1, 1) AS ix) tuple
 223       ORDER BY ix) inn;
 224 $$ LANGUAGE SQL IMMUTABLE STRICT;
 225
 226 create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
 227 declare
 228         adj integer[];
 229 begin
 230         if array_length(hex,1) = 2 then
 231                 hex = hex.ijk(hex);
 232         end if;
 233
 234         dir := dir % 6;
 235         if dir = 2 then
 236                 hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
 237         elsif dir = 1 then
 238                 hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
 239         elsif dir = 0 then
 240                 hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
 241         elsif dir = 5 then
 242                 hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
 243         elsif dir = 4 then
 244                 hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
 245         elsif dir = 3 then
 246                 hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
 247         else
 248                 -- error
 249         end if;
 250         return hex;
 251 end;
 252 $$ immutable strict;
 253
 254 create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
 255         select hex.adjhex(dir, hex.ijk(hex));
 256 $$ immutable strict;
 257
 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);
 261 $$ immutable strict;
 262
 263 create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
 264 begin
 265         -- convert to xy
 266         if array_length(hex,1) = 3 then
 267                 hex = hex.xy(hex);
 268         end if;
 269
 270         if hex[1] < 0 or hex[2] < 0 then
 271                 hex[1] = hex[1] * 2;
 272                 hex[2] = hex[2] * 2;
 273                 if hex[1] < 0 then
 274                         hex[1] = -hex[1] - 1;
 275                 end if;
 276                 if hex[2] < 0 then
 277                         hex[2] = -hex[2] - 1;
 278                 end if;
 279                 return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 280         end if;
 281
 282         return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
 283 end;
 284 $$ immutable strict;
 285
 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;
 288 $$ immutable strict;
 289
 290 create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
 291         select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
 292 $$ immutable strict;
 293
 294 create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
 295 declare
 296         c integer[]; -- cantor coordinates
 297         base integer[];
 298 begin
 299         if range <= 0 then
 300                 return;
 301         end if;
 302         if include then
 303                 return next hex;
 304         end if;
 305
 306         -- we could in theory eliminate the inner loop and
 307         -- do a return query
 308         base = hex.ijk(hex);
 309         for r in 1 .. range loop
 310                 c = base;
 311                 c[1] := c[1] + r;
 312                 c[3] := -c[1] - c[2]; -- keep the invariant
 313                 hex := hex.cantor(c);
 314
 315                 for h in 0 .. r * 6 - 1 loop
 316                         return next hex;
 317                         hex = hex.adjacent(h/r+2,hex);
 318                 end loop;
 319         end loop;
 320         return;
 321 end;
 322 $$ immutable strict;
 323
 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;
 326 $$ immutable strict;
 327
 328 create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
 329 declare
 330         c integer[]; -- cantor coordinates
 331 begin
 332         if range = 0 then
 333                 return next hex;
 334         end if;
 335         if range <= 0 then
 336                 return;
 337         end if;
 338
 339         c = ijk(hex);
 340         c[1] := c[1] + range;
 341         c[3] := -c[1] - c[2]; -- keep the invariant
 342         hex := hex.cantor(c);
 343
 344         for q in 0 .. range * 6 - 1 loop
 345                 return next hex;
 346                 hex = hex.adjacent(q/range+2,hex);
 347         end loop;
 348         return;
 349 end;
 350 $$ immutable strict;
 351
 352 select
 353 xy(ARRAY[1,2]),
 354 ijk(xy(ARRAY[1,2])),
 355 xy(ijk(xy(ARRAY[1,2])))
 356 ;
 357
 358 select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
 359 select
 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;
 364
 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;