Tile Serving with Dynamic Geometry
In our previous posting on tile serving, we showed off how pg_tileserv can use database functions to generate tiles by querying existing tables with user parameters.
We can also use functions to build geometry on the fly without input from tables. For example, hexagons!
Hexagons are a favourite input for visualizations, because they have a nice symmetric shape and provide equal areas for summarization.
Hexagonal Grid
A filling of the plane with hexagons is a hexagonal grid. Curiously, it's possible to address any hexagon in the grid with a two-dimensional coordinate, just like with a square grid, the only difference being that every second column of the hexgrid is offset one-half a hexagon vertically.
Generating a Hexagon
To build a hexagonal tiling we need a function that produces an arbitrary hexagon in a hexagon grid, like this.
CREATE OR REPLACE
FUNCTION
hexagon(i integer, j integer, edge float8)
RETURNS geometry
AS $$
DECLARE
h float8 := edge*cos(pi()/6.0);
cx float8 := 1.5*i*edge;
cy float8 := h*(2*j+abs(i%2));
BEGIN
RETURN ST_MakePolygon(ST_MakeLine(ARRAY[
ST_MakePoint(cx - 1.0*edge, cy + 0),
ST_MakePoint(cx - 0.5*edge, cy + -1*h),
ST_MakePoint(cx + 0.5*edge, cy + -1*h),
ST_MakePoint(cx + 1.0*edge, cy + 0),
ST_MakePoint(cx + 0.5*edge, cy + h),
ST_MakePoint(cx - 0.5*edge, cy + h),
ST_MakePoint(cx - 1.0*edge, cy + 0)
]));
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE
STRICT
PARALLEL SAFE;
Any unique edge length produces a unique gridding of the plane, so to build a hexagon, we need both the hexagon coordinates, and the edge length that defines the gridding.
Filling a Tile with Hexagons
Each tile will need to know what hexagons to generate, from the overall hexagon grid, so we also need a function that can take a tile bounds and figure out which hexagons fall inside the bounds.
CREATE OR REPLACE
FUNCTION hexagoncoordinates(bounds geometry, edge float8,
OUT i integer, OUT j integer)
RETURNS SETOF record
AS $$
DECLARE
h float8 := edge*cos(pi()/6);
mini integer := floor(st_xmin(bounds) / (1.5*edge));
minj integer := floor(st_ymin(bounds) / (2*h));
maxi integer := ceil(st_xmax(bounds) / (1.5*edge));
maxj integer := ceil(st_ymax(bounds) / (2*h));
BEGIN
FOR i, j IN
SELECT a, b
FROM generate_series(mini, maxi) a,
generate_series(minj, maxj) b
LOOP
RETURN NEXT;
END LOOP;
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE
STRICT
PARALLEL SAFE;
This set-returning function produces all the hexagon coordinates that interact with a bounds in the hex plane.
So, for example, given the bounds of a tile, it generates a set of records with hex coordinates:
SELECT * FROM hexagoncoordinates(ST_TileEnvelope(15, 1, 1), 1000.0);
i | j
--------+-------
-13358 | 11567
-13358 | 11568
-13357 | 11567
-13357 | 11568
-13356 | 11567
-13356 | 11568
Generating Hexagons for Any Tile
Remember, for pg_tileserv
to publish a function as a tile set, the function needs to take in the tile coordinates (x, y and zoom) as well as any other user parameters (in our case, a "step") parameter.
The "step" parameter just controls the hexagon edge size, as a function of the tile width. A "step" of 1 implies an edge size 1/2 of the tile size; a "step" of 2 implies an edge size 1/4 of the tile size, and so on.
CREATE OR REPLACE
FUNCTION tilehexagons(z integer, x integer, y integer, step integer,
OUT geom geometry(Polygon, 3857), OUT i integer, OUT j integer)
RETURNS SETOF record
AS $$
DECLARE
bounds geometry;
maxbounds geometry := ST_TileEnvelope(0, 0, 0);
edge float8;
BEGIN
bounds := ST_TileEnvelope(z, x, y);
edge := (ST_XMax(bounds) - ST_XMin(bounds)) / pow(2, step);
FOR geom, i, j IN
SELECT ST_SetSRID(hexagon(h.i, h.j, edge), 3857), h.i, h.j
FROM hexagoncoordinates(bounds, edge) h
LOOP
IF maxbounds ~ geom AND bounds && geom THEN
RETURN NEXT;
END IF;
END LOOP;
END;
$$
LANGUAGE 'plpgsql'
IMMUTABLE
STRICT
PARALLEL SAFE;
Generating the Tile
Now that we can generate a set of hexagons for any tile coordinate, we just need a function that wraps up those hexagons in an MVT format for output.
CREATE OR REPLACE
FUNCTION public.hexagons(z integer, x integer, y integer, step integer default 4)
RETURNS bytea
AS $$
WITH
bounds AS (
-- Convert tile coordinates to web mercator tile bounds
SELECT ST_TileEnvelope(z, x, y) AS geom
),
rows AS (
-- All the hexes that interact with this tile
SELECT h.i, h.j, h.geom
FROM TileHexagons(z, x, y, step) h
),
mvt AS (
-- Usual tile processing, ST_AsMVTGeom simplifies, quantizes,
-- and clips to tile boundary
SELECT ST_AsMVTGeom(rows.geom, bounds.geom) AS geom,
rows.i, rows.j
FROM rows, bounds
)
-- Generate MVT encoding of final input record
SELECT ST_AsMVT(mvt, 'public.hexagons') FROM mvt
$$
LANGUAGE 'sql'
STABLE
STRICT
PARALLEL SAFE;
COMMENT ON FUNCTION public.hexagons IS 'Hex coverage dynamically generated. Step parameter determines how approximately many hexes (2^step) to generate per tile.';
The result is an infinitely zoomable hexagon grid that is generated on-the-fly.
Related Articles
- PostGIS Day 2024 Summary
8 min read
- Crunchy Data Warehouse: Postgres with Iceberg for High Performance Analytics
8 min read
- Loading the World! OpenStreetMap Import In Under 4 Hours
6 min read
- Easy Totals and Subtotals in Postgres with Rollup and Cube
5 min read
- A change to ResultRelInfo - A Near Miss with Postgres 17.1
8 min read