db/functions/functions.sql

   1 --------------------------------------------------------------------------------
   2 -- SQL versions of the C database functions.
   3 --
   4 -- Pure pl/pgsql versions are *slower* than the C versions, and not recommended
   5 -- for production use. However, they are significantly easier to install, and
   6 -- require fewer dependencies.
   7 --------------------------------------------------------------------------------
   8
   9 -- tile_for_point function returns a Morton-encoded integer representing a z16
  10 -- tile which contains the given (scaled_lon, scaled_lat) coordinate. Note that
  11 -- these are passed into the function as (lat, lon) and should be scaled by
  12 -- 10^7.
  13 --
  14 -- The Morton encoding packs two dimensions down to one with fairly good
  15 -- spatial locality, and can be used to index points without the need for a
  16 -- proper 2D index.
  17 CREATE OR REPLACE FUNCTION tile_for_point(scaled_lat int4, scaled_lon int4)
  18   RETURNS int8
  19   AS $$
  20 DECLARE
  21   x int8; -- quantized x from lon,
  22   y int8; -- quantized y from lat,
  23 BEGIN
  24   x := round(((scaled_lon / 10000000.0) + 180.0) * 65535.0 / 360.0);
  25   y := round(((scaled_lat / 10000000.0) +  90.0) * 65535.0 / 180.0);
  26
  27   -- these bit-masks are special numbers used in the bit interleaving algorithm.
  28   -- see https://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN
  29   -- for the original algorithm and more details.
  30   x := (x | (x << 8)) &   16711935; -- 0x00FF00FF
  31   x := (x | (x << 4)) &  252645135; -- 0x0F0F0F0F
  32   x := (x | (x << 2)) &  858993459; -- 0x33333333
  33   x := (x | (x << 1)) & 1431655765; -- 0x55555555
  34
  35   y := (y | (y << 8)) &   16711935; -- 0x00FF00FF
  36   y := (y | (y << 4)) &  252645135; -- 0x0F0F0F0F
  37   y := (y | (y << 2)) &  858993459; -- 0x33333333
  38   y := (y | (y << 1)) & 1431655765; -- 0x55555555
  39
  40   RETURN (x << 1) | y;
  41 END;
  42 $$ LANGUAGE plpgsql IMMUTABLE;
  43
  44
  45 -- maptile_for_point returns an integer representing the tile at the given zoom
  46 -- which contains the point (scaled_lon, scaled_lat). Note that the arguments
  47 -- are in the order (lat, lon), and should be scaled by 10^7.
  48 --
  49 -- The maptile_for_point function is used only for grouping the results of the
  50 -- (deprecated?) /changes API call. Please don't use it for anything else, as
  51 -- it might go away in the future.
  52 CREATE OR REPLACE FUNCTION maptile_for_point(scaled_lat int8, scaled_lon int8, zoom int4)
  53   RETURNS int4
  54   AS $$
  55 DECLARE
  56   lat CONSTANT DOUBLE PRECISION := scaled_lat / 10000000.0;
  57   lon CONSTANT DOUBLE PRECISION := scaled_lon / 10000000.0;
  58   zscale CONSTANT DOUBLE PRECISION := 2.0 ^ zoom;
  59   pi CONSTANT DOUBLE PRECISION := 3.141592653589793;
  60   r_per_d CONSTANT DOUBLE PRECISION := pi / 180.0;
  61   x int4;
  62   y int4;
  63 BEGIN
  64   -- straight port of the C code. see db/functions/maptile.c
  65   x := floor((lon + 180.0) * zscale / 360.0);
  66   y := floor((1.0 - ln(tan(lat * r_per_d) + 1.0 / cos(lat * r_per_d)) / pi) * zscale / 2.0);
  67
  68   RETURN (x << zoom) | y;
  69 END;
  70 $$ LANGUAGE plpgsql IMMUTABLE;