Tiles are not only pivotal for visualizing map data in browsers, but they can also revolutionize the way we store geospatial data. This article explores the use of a ’tile-based GIS’ for storing data like population distribution, utilizing PostGIS to efficiently export this data into vector tiles (MVT) for web-based dissemination."
Certain datasets, such as global population distribution or OpenStreetMap buildings’ completeness assessments, are naturally suited to a grid or ’tile-based GIS’ format. This format, relying on either coordinate system grids (like a 30 arc-seconds grid) or tile concepts, can be adeptly managed using the quadtree structure with quadkeys.
This blog post focuses on addressing tile-based data using quadkeys and demonstrates how to serve this data efficiently using PostgreSQL’s PostGIS extension, specifically its ST_asMVT function.
Understanding Tiles
Tiles refer to partitioning a two-dimensional space, such as a projected Earth, by recursively subdividing it into quadrants. Each subdivision step, or zoom level, further divides each quadrant into smaller ones. Tiles are typically addressed using simple coordinates (x, y) along with the zoom level (z).
Quadkeys
Quadkeys, based on the quadtree structure, offer a convenient method for addressing tiles. Each quadkey digit represents a zoom level, making them intuitive for both reading and applying aggregation logic.
PostgreSQL/PostGIS Functions
The following PostgreSQL/PostGIS functions facilitate the conversion between XYZ tile coordinates and quadkeys:
---- CONVERT XYZ TO QUADKEY
Especially for the visualization they became
CREATE OR REPLACE FUNCTION public.convert_xyz_to_quadkey(x INTEGER, y INTEGER, z INTEGER)
RETURNS TEXT AS
$$
DECLARE
quadkey TEXT;
quadkey_digit INTEGER;
rx INTEGER;
ry INTEGER;
BEGIN
FOR i IN 1..z LOOP
rx = x % 2;
x = x / 2;
ry = y % 2;
y = y / 2;
IF rx = 0 AND ry = 0 THEN
quadkey_digit = 0;
ELSIF rx = 1 AND ry = 0 THEN
quadkey_digit = 1;
ELSIF rx = 0 AND ry = 1 THEN
quadkey_digit = 2;
ELSIF rx = 1 AND ry = 1 THEN
quadkey_digit = 3;
ELSE
RAISE EXCEPTION 'Failure converting xyz to Quadkey';
END IF;
quadkey = CONCAT(CAST(quadkey_digit AS CHAR), quadkey);
END LOOP;
RETURN quadkey;
END;
$$
LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE;
COMMENT ON FUNCTION public.convert_xyz_to_quadkey IS 'Converts a tile identifier from the format x,y,z to a Quadkey.';
---- CONVERT QUADKEY TO XYZ
CREATE OR REPLACE FUNCTION public.convert_quadkey_to_xyz(quadkey TEXT)
RETURNS SETOF INTEGER AS
$$
DECLARE
x INTEGER := 0;
y INTEGER := 0;
z INTEGER := 0;
quadkey_digit TEXT;
BEGIN
FOR i IN 1..LENGTH($1) LOOP
z = z + 1;
quadkey_digit = SUBSTRING(quadkey, i, 1);
IF quadkey_digit = '0' THEN
x = x * 2; y = y * 2;
ELSIF quadkey_digit = '1' THEN
x = x * 2 + 1; y = y * 2;
ELSIF quadkey_digit = '2' THEN
x = x * 2; y = y * 2 + 1;
ELSIF quadkey_digit = '3' THEN
x = x * 2 + 1; y = y * 2 + 1;
ELSE
RAISE EXCEPTION 'Failure: a Quadkey can only contain numbers from `0` to `3`.';
END IF;
END LOOP;
RETURN QUERY VALUES (x), (y), (z);
END
$$
LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE;
COMMENT ON FUNCTION public.convert_quadkey_to_xyz IS 'Converts a tile identifier from the format Quadkey to x,y,z.';
Serving Tiles with PostGIS
To serve this data as vector tiles over the internet, we use the ST_AsMVT function in PostGIS. This function can efficiently transform spatial data stored in PostGIS into vector tiles that can be rendered in a web map.
---- Generates tile geometry from Quadkey
CREATE OR REPLACE FUNCTION public.generate_tile_geometry_from_quadkey(quadkey TEXT)
RETURNS GEOMETRY -- Using Web Mercator, EPSG 3857
AS $$
DECLARE
earth_rad FLOAT; x FLOAT; y FLOAT; off_set FLOAT; q TEXT;
BEGIN
SELECT 6378137.0 INTO earth_rad; -- EPSG 3857 definition
SELECT 0 INTO x;
SELECT 0 INTO y;
SELECT earth_rad * pi() INTO off_set; -- Tile centroid
FOR i IN 1..LENGTH($1) LOOP
q = SUBSTRING(quadkey, i, 1);
off_set = off_set * .5;
IF q = '0' THEN
x = x - off_set; y = y + off_set;
ELSIF q = '1' THEN
x = x + off_set; y = y + off_set;
ELSIF q = '2' THEN
x = x - off_set; y = y - off_set;
ELSIF q = '3' THEN
x = x + off_set; y = y - off_set;
ELSE RAISE EXCEPTION 'Failure: a Quadkey can only contain numbers from `0` to `3`.';
END IF;
END LOOP;
RETURN ST_MakeEnvelope(x - off_set, y - off_set, x + off_set, y + off_set, 3857);
END
$$
LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE;
---- Generates vector tiles containing their own geometry
CREATE OR REPLACE FUNCTION public.generate_vector_tile_with_geometry(z INTEGER, x INTEGER, y INTEGER)
RETURNS BYTEA AS
$$
DECLARE
vector_tile BYTEA DEFAULT '';
BEGIN
WITH mvtdata AS (
SELECT ST_MakeEnvelope(0, 0, 4096, 4096) AS geom,
z as zoom
)
SELECT ST_AsMVT(mvtdata.*, 'tiles' 4096, 'geom') FROM mvtdata INTO vector_tile;
RETURN vector_tile;
END;
$$
LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE;
---- Generates vector tiles containing geometries of tiles with a specified zoom level and allows for vizualization on lower zoom levels
CREATE OR REPLACE FUNCTION public.generate_vector_tile_containing_geometries_of_sub_tiles(z INTEGER, x INTEGER, y INTEGER, target_zoom_level INTEGER)
RETURNS BYTEA AS
$$
DECLARE
vector_tile BYTEA DEFAULT '';
inside_tiles INTEGER; -- Number of tiles at zoom level `target_zoom_level` along one axis of a tile with zoom level `z`.
inside_tile_pixels INTEGER;
BEGIN
IF z > target_zoom_level THEN
RAISE NOTICE 'Vector tiles cannot be be created for zoom levels higher than the indicated in `target_zoom_level`: %.', target_zoom_level;
ELSE
inside_tiles = 2 ^ (target_zoom_level - z);
inside_tile_pixels = 4096 / inside_tiles;
WITH mvtdata AS (
SELECT ST_MakeEnvelope(
inside_tile_pixels * tile_x, inside_tile_pixels * tile_y,
inside_tile_pixels * (tile_x + 1), inside_tile_pixels * (tile_y + 1)
) AS geom
FROM
generate_series(0, inside_tiles) AS tile_x,
generate_series(0, inside_tiles) AS tile_y
)
SELECT ST_AsMVT(mvtdata.*, 'tiles', 4096, 'geom') FROM mvtdata INTO vector_tile;
END IF;
RETURN vector_tile;
END;
$$
LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE;